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CHAPTER  I 
INTRODUCTION 

1.1  Motivation  and  General  Purpose 

In  the  determination  of  the  amount  of  reflected  electro- 
magnetic energy  from  a  surface,  called  the  radar  backscattering, 2 
one  of  the  major  methods  of  modeling  the  reflecting  surface  is 
to  assume  the  surface  height  above  some  mean  surface  can  be  de- 
scribed by  a  random  process.9   By  the  tangent  plane  approximation 
which  assumes  the  reflected  body  is  isotropic  and  homogeneous, 
the  radar  backscatter  can  be  obtained  in  an  integral  form. 
Essentially  the  integrand  of  this  integral  is  a  new  random  pro- 
cess which  is  a  function  of  the  surface  height  random  process 
and  its  derivatives.   The  form  of  this  function  changes  as  the 
shape  of  the  average  surface  changes,  i.e.,  a  sphere  or  a  plane. 
It  is  desired  to  find  the  probability  density  of  the  integrand 
and  the  integral  for  each  case. 

Upon  the  realization  of  the  magnitude  of  this  problem,  a 
sub-problem  is  decided  upon  which  could  be  utilized  as  an  interim 
step  regardless  of  the  form  of  the  average  surface.   The  inte- 
grands in  most  cases  can  be  expressed  in  terms  of  a  power  series 
of  cos  0-^  times  an  exponential  term  whose  exponent  is  a  function 
of  only  the  surface  height,  where  Q±    is  the  angle  between  the 
surface  normal  and  the  direction  to  the  receiver.   Therefore  as 
a  first  step,  it  is  decided  to  determine  the  probability  density 
of  cos  02  and  then  this  density  can  later  be  applied  to  a  par- 
ticular surface. 


1.2   The  Problem 

For  the  problem  of  radar  backscatterlng  from  a  rough  spheri- 
cal  surface,   normal  spherical  coordinates  (a   ,  aQ,  aw-)  are 
defined.   The  variation  in  the  radius  of  the  rough  sphere  from 
the  average  at  a  point  (0,  0)    is  denoted  by  H(0,  0)  .   The  random 
surface  is  assumed  continuous  in  the  mean  and ' diff erentiable  over 
a  finite  region.   The  radius  vector  from  the  origin  to  the  sur- 
face  point  is  rQ.   (rQ  =  rQ  ar  ).   Then  the  equation  of  the 
surface  can  be  written  as 


^  =  rQ  -  [(a  +  H(e,j2T)]  =  0 


(1.1.1) 


where  a  is  an  average  radius  of  the  rough  surface.   By  normaliz- 
ing the  gradient  of  the  equation  ( 1.1.1),  the  unit  outward  sur- 


face normal  vector  a   is  found. 
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Performing  the  indicated  operations  of  (1.1.2)  on  ^as  defined 
in  (1.1.1),  the  above  equation  yields 
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(1.1. k) 


Under  the  assumption  that  H(0,  0)    is  a  normal  variable,  its 
partial  derivatives  with  respect  to  its  variables  are  also 
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normal.   '    Therefore  two  normal  variables  are  defined  as 

1   d 

x  =  —  —  H(e,  $) 
r0  ae 

1        d 

Y  =  _  H(e,  jzf) 

r0  sin  9  d® 

aR  denotes  the  direction  of  the  receiver  from  the  point  on 
the  rough  surface  given  by  rQ,  6,  J2f.   The  angle  between  a   and 
an  is  defined  as  6  .   The  expression  of  cos  6,  then  is  found. 

cos  6  +  Y  sin  6 

cos  ei  =  IT        9   o  (i.i.5) 

-Y  1  +  X2  +  Y2 

Since  cos  6-^  involves  the  function  fl  +  X2  +  Y2,  the 
problem  of  finding  the  density  function  of  cos  6-,  starts  from 
finding  the  densities  of  71  +  X2  +  y2  and  l/Vl  +  X2  +  Y2  with 
a  desire  to  find  the  former  density  by  the  information  contained 
in  the  latter  two. 

Hopefully,  the  probability  density  of  cos  B-.    can  be  deter- 
mined analytically  and  the  variations  of  this  density  with  each 
different  standard  deviation  a~  ,    correlation  coefficient  f,    and 
angle  6  can  be  studied.   If  it  fails,  a  valid  approximation  is 
desired.   The  Monte  Carlo  method  is  thus  introduced  to  solve  the 
aforementioned  problem  numerically  due  to  the  failure  of  finding 
the  density  function  of  cos  6^  analytically. 

1.3   Summary  of  Chapter  Development 

Chapter  I  contains  the  statement  of  the  problem.   Chapter  II 


presents  two  general  analytical  methods  for  finding  the  prob- 
ability density  function.   The  explicit  densities  of 
Vl  +  X2  +  Y2  and  l/Vl  +  X2  +  Y2  are  obtained  by  using  these 
two  methods.   Appendix  A  presents  the  derivation  of  the  density 
function  of  cos  9-j_  by  analytical  methods.   However,  the  problem 
is  not  solved  by  these  methods. 

Chapter  III  states  the  Monte  Carlo  method  used  to  obtain 
the  desired  probability  density.   First,  the  method  of  generating 
random  numbers  is  discussed.   Chi-square  test  is  mentioned  here 
to  check  the  goodness  of  fit  of  the  generated  random  numbers. 
Two  curve  fitting  techniques  are  developed  for  finding  the 
explicit  expression  by  the  given  data  from  the  Monte  Carlo  method. 
More  details  of  the  curve  fitting  algorithm  are  given  in 
Appendix  B. 

Chapter  IV  presents  the  results  of  densities  by  using  the 
methods  of  Chapter  III.   A  comparison  is  made  to  the  density 
function  of  i4  +  X2  +  Y2  and  l/Vl  +  X2  +  Y2  with  the  results 
by  the  analytical  methods.   The  densities  of  cos  Bi   for  certain 
parameters  are  also  obtained.   Appendix  C  is  the  computer  pro- 
grams used  in  this  chapter.   The  various  densities  of  cos  6-, 
corresponding  to  different  parameters  are  discussed  and  further 
researches  are  recommended  in  Chapter  V. 


CHAPTER  II 

THE  ANALYTICAL  METHODS  OP  FINDING  THE  DENSITY 
FUNCTION  OF  A  FUNCTION  OF  TWO  OR 
MORE  RANDOM  VARIABLES 

2.1   Introduction 

The  difficulty  of  finding  the  probability  density  of  a 
function  of  random  variables  depends  on  the  number  of  random 
variables  and  the  complexity  of  the  functional  relationships 
involved.   If  these  random  variables  are  independent,  it  makes 
the  problem  easier  than  if  they  are  dependent.   The  densities 
can  be  quickly  calculated  for  linear  functions  of  random  vari- 
ables.  In  most  physical  problems,  either  the  random  variables 
are  not  independent  or  the  function  is  nonlinear.   When  either 
or  both  of  these  cases  occur  the  algebraic  solutions  require 
very  complicated  integrations  which  are  not  always  solvable  in 
closed  form  or  expressed  as  a  well-known  series. 

In  this  chapter  two  general  analytical  methods  of  finding 
the  density  function  will  be  discussed.   The  algebraic  solutions 
for  the  density  functions  of  Vl  +  X2  +  Y2  and  l/Vl  +  X2  +  Y2 
are  obtained  by  these  methods.   These  explicit  densities  will  be 
compared  with  the  densities  obtained  by  the  Monte  Carlo  method 
in  order  to  estimate  the  accuracy  and  validity  of  the  Monte 
Carlo  method. 


2.2   The  Density  Function  of  a  Function  of  Two  or 

More  Random  Variables 

Under  certain  conditions,  ^>   the  function  of  random  vari- 
ables can  be  considered  as  another  random  variable.   There  are 
many  methods  of  finding  the  density  function  of  this  new  random 
variable.   Only  two  different  methods  are  discussed  here. 

The  first  method  consists  of  finding  the  distribution  func- 
tion Fz(z)  first,  then  by  taking  derivative  with  respect  to  z 
the  probability  density  function  is  obtained. 

fz  (z)  =  —  Fz  (z)  (2.2.1) 

In  order  to  determine  Fz(z)  for  a  given  z,  the  probability  of 
the  event  \Z   ^  zj  must  be  found.   For  the  function  of  two  random 
variables  case,  Dz  is  denoted  by  the  region  of  XY  plane  such 
that  g(x,y)  =  z.   Then  (z  ^   z)  =  ((X,Y)  £  Dzj  .   Hence  it  suf- 
fices to  find  the  probability  mass  in  the  region  Dz.   This  mass 
is  given  by  the  integration 

Fz(z)  =  P(Z  ^  z)  =  p  (X,  Y)  £   D7 


-if 


^XY  (x>  ^dx  ^  (2.2.2) 

Dz' 

where  fXy(x,  y)  is  the  original  joint  density  function  of 
random  variables  X  and  Y. 

Similarly  for  a  function  of  n  random  variables,  that  is 
z  =  g(x1,  x2,  .  .  .,  xn) ,  Dz  is  denoted  by  the  region  of  n- 
dimensional  space  such  that  g(xx,  x2,  .  .  . ,  xn)  =  z  and 


{z  ^   zj    ={(X1,    X2,    .    .    .,    Xn)   £   Dzj  .      Then 

(z)    =  P(Z£    z)    =  p{(X1,    X2,    .    .  .,    Xn)    €   Dz] 

=    J      ...|fYY  (xn,  x2,    ...,    xn)dXldx2, 

>DZ        /      A1A2    * ' '    An 


FZ> 


dxn 


(2.2.3) 


where   f x  x      x   (x1,  x2,  ..-,  *n)  is  the  original  joint 

density  function  of  random  variables  X-.,  X2,  .  .  .,  X  . 

The  second  method  consists  of  introducing  an  auxiliary  vari- 
able and  finding  the  joint  density  function  of  the  desired 
random  variable  and  the  auxiliary  variable.   Then  by  integrating 
the  joint  density  function  over  the  entire  range  of  the  auxiliary 
variable,  the  marginal  density  function  of  the  random  variable 
is  obtained. 

For  the  case  of  the  function  of  two  random  variables,  the 
auxiliary  variable  ¥  is  assumed  to  be  W  =  g2(X,  Y) .   With 
Z  =  g-i(X,  Y)  the  joint  density  is 

x    fXY(xl>yi)          fXY(xn>Yn)         ,     .  , 
f™  (z,w)  =  -+  .  .  .  + (2.2.I4.) 


ZW 


J(xi,yi)|  lJ(xn*yn)l 


where  (x1,y1),  (x2,y'2)  .  .  .  (xn,yn)  are  all  the  real  solutions 
of  the  equation,  g-,(x,y)  =  z  with  g2(x,y)  =  w,  and  J(x,y)  is 
the  Jacobian  of  the  transformation.   The  unknown  density  func- 
tion of  Z  is 

(°° 

fz(z)  =      f^z.wjdw  (2.2.5) 

The  second  method  can  be  extended  to  a  function  of  n  random 
variables.   Given  Z-,  =  g-,(X-,,  X2,  ...,  X  ),  then  n  -  1  auxiliary 
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variables  as 

z2  =  g2(x1,x2,...,xn),  z3  =  g3(x1,x2,...,xn),  .  .  ., 

Zn  =  Sn(xi>X2>- • ->Xn) 

are  assumed.   The  choice  of  the  auxiliary  variables  depends  on 

which  assumption  can  simplify  the  problem  of  finding  the  joint 

density.   By  solving  these  n-system  equations,  k  sets  of  real 

solution  (x    x    x    ...,  x  .)  are  obtained,  where 
-Lx   d.x        ix  ni 

1  =  1,  2,  3,  .  .  .,  k.   Therefore  the  joint  density  function  is 

f         (*      „        \         £-       fXiX2...Xn  (xli>x2i>- • -'Xni) 

ZlZ2'"zn  X   2""'Zn  =  ?    |~^ 

i-1   |J(xli,x2i,. ..,  ...,..., xni) 

(2.2.6) 
where  J(x±t    xg,  ...,  xn)  is  the  Jacobian  of  the  transformation 
of  these  n-system  equations.   Then 

fzi(Zl)  =/-oo"*  |.00fzlZ2...Zn(Z^Z2---Zn)dz2dz3  '  '  '  dzn 

.  (2.2.7) 

The  above  two  methods  are  commonly  used  to  find  the  density 
function  of  a  function  of  random  variables.   They  will  be  used 
in  the  next  section. 

2.3   The  Density  Functions  of  J 1  +  X2  +  Y2  and  l/Vl  +  X2  +  Y2 

For  finding  the  density  functions  of  -f/ 1  +  X2   +  Y2   and 
l/yl  +  X2  +  y2,  the  polar  coordinates  are  used  here  due  to  their 
simplicity  in  the  operation.   Assuming  Z  =  X2  +  Y2,  the  density 
function  of  Z  is  obtained  by  the  first  method  of  the  last  sec- 
tion.  After  that,  a  new  random  variable  W  is  defined  as  Z  +  1. 


Its  density  is  found  merely  by  shifting  the  variable  of  the 
density  function  of  Z  to  If  -  1,   The  square  root  of  W  and  one 
over  this  square  root  are  other  random  variables.   These  two 
densities  are  obtained  by  the  second  method  of  the  last  section 
based  on  the  density  of  W. 

X  and  Y  are  two  random  variables  with  zero  means  and  equal 
standard  deviation  <T   and  correlation  coefficient  /° .      The  joint 
probability  density  function  of  X  and  Y  is 

?T£^>?)    =  ~ 2  j  g  exp  ( _       (2.3.1) 

By   starting  with  Z  =  X2   +  Y2,    then  Dz   is   a    circle  with  radius 
of  ATz.      In   polar  form,    r  =   YT,    x  =   r  cos   0,    y  =  r   sin   6,    then 

PZ(z)    =  (     fexpf _    .    (x2+y2-2/4cy)]dxdy 

1  f^  f         -r2  "\    ,2%         (p{sln  29)r2) 

=  2„r2  fT^    i        '   eXP(2^2(lV2)j  I      eXP(2CT-2(1.^2)  )  dMr 

The  second  part  of  the  above  integration  is  the  modified  Bessel 
function  of  first  kind  of  order  zero1  denoted  as  I  ;  then 


M^P*    }0 


yr      (-    -r2     )    s    pr 


2 


By  taking  the  derivative  of  Pz(z)  with  respect  to  z  and  using 
the  generalized  Leibnitz  formula,^  the  density  function  of  Z 
is  obtained. 
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-z 


Msb)  = 


exp- 


2cr2/]772    [20-2(1-/32) 


*o/* 


2cr2(i-/^2) 


u(z) 
(2.3.2) 


where  u(z)  is  the  unit  step  function. 

Defining  W  =  Z  +  1,  with  fw(w)  =  fz(w  -  1),  then 


%(w) 


2cr2Vir72 


/-   -(w-l) 
exp-f 1  I 


/>(w-l) 


20-2(1-  /)2)J   °)  2CT2(l-/02) 


u(w-l) 


(2.3.3) 


For 


S  =7W  =T1  +  X2  +  Y2  ,  with  fs(s)  =  2s  f¥(s2),  then 


2s 
fsU)  =  20-2^1772 


exp 


(     -(az-l)   -)    ,-     (s^-1)   ) 
2Cr2(i-/>2)J   °^  20-2(1-  /72,J 


(2.3.4) 


For  V  =  - 


S   4/l+X2+Y2^ 


1      1 
with  fy(v)  =  —  fs( ),  then 


fv(v)  = 


,    -d-v2)     -J     •    d-v^)r    ) 

eXPl2cr2(1_/)2)v2J  I°l2cr2(i-/^2)v2j 


[u(v)  -  u(v  -  1)] 


(l-v2^ 

(2.3.5) 


Equations  (2. 3. 4)  and  (2.3.5)  are  the  solutions  needed  to  com- 
pare with  the  results  of  the  Monte  Carlo  method  in  Chapter  IV. 
They  are  also  plotted  for  <T  -   1  and  r  -  ®-k->    0.9  as  Fig.  6 
through  Fig.  9  in  Chapter  IV  for  a  comparison. 

For  the  limiting  case,  when  /^  =  1,  equations  ( 2.3.1]-)  and 
(2.3.5)  yield 


fs(s)  = 


<7-y-n(s2-l)      I  \±CT 2 


-(s2-l) 
exp  } fu(s  -  1) 


(2.3.6) 
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1                           /     (1-vO)    r  -, 

fv(v)    =  exp —      fu(v)    -  u(v-l)J         (2.3.7) 

When     p=  0,    they  yield 

s  f     (s2   -    1)"\ 

fS(s)    =  — o   exP1~  5— /  u(s   ~   1}  (2.3.8) 

(T^  (  2CT2     J 

1  (     (l-v2))r 

fV(v)   =  — T"J  exP  1 TT  f  I  u^v)    "  u<v  "   Dj  (2.3.9) 

cr2v3  [     2<Tlvl)  L  J 

The  graphs  of  equations  (2.3.6)  -  (2.3.9)  are  shown  in  Fig.  1 
and  Fig.  2  for  <7~  =  1  case. 

All  the  results  of  the  above  density  functions  are  verified 
by  integrating  over  the  entire  range  of  the  variable  and  obtain- 
ing the  number  one. 

The  derivation  of  the  density  function  of  cos  6n  is  shown 
in  Appendix  A.   With  the  methods  used  here  and  the  method  by 
Fourier  transformation,  all  involve  a  difficult  integration. 
It  cannot  end  with  a  closed  form  or  a  well-known  series  expres- 
sion.  Therefore  the  Monte  Carlo  method  is  used  in  Chapter  IV 
which  gives  an  approximate  result  with  certain  fixed  parameters. 
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Pig.  1.   Density  function  of  S  with  g-  =  1  by  the 

Analytical  method  (S  =  >j/l  +  X2  +  Y2). 
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Pig.  2.   Density  function  of  V  with  <7~  =  1  by  the 
analytical  method  (V  =  l/Yl  +   X2  +  Y2)  . 


CHAPTER  III 
THE  MONTE  CARLO  METHOD  AND  CURVE  PITTING  TECHNIQUE 

3.1   Introduction 

As  previously  indicated,  analytical  methods  of  obtaining 
probability  densities  for  functions  of  random  variables  do  not 
always  yield  closed  form  solutions.   In  this  chapter,  the  Monte 
Carlo  method  is  introduced.   When  it  is  combined  with  curve 
fitting  techniques,  an  approximate  closed  form  solution  can 
be  obtained. 

The  important  job  of  the  Monte  Carlo  method  is  to  generate 
a  set  of  random  numbers  with  a  certain  distribution.   How  to 
generate  random  numbers  of  uniform  distribution  and  how  to  change 
them  to  another  distribution  are  discussed  in  section  3.2.   This 
job  will  be  done  by  the  computer,  IBM  36o/50.   With  a  limited 
sampled  size  and  computer  round-off  errors,  the  accuracy  of  the 
Monte  Carlo  method  has  to  be  considered.   It  is  discussed  in 
section  3-3>- 

Finally,  two  curve  fitting  techniques  which  will  be  used 
to  find  the  density  function  explicitly  with  the  data  from  the 
Monte  Carlo  method  are  developed  in  section  3.1].. 

3.2  Generation  of  Random  Numbers 

The  Monte  Carlo  method10'12'1^  is  a  kind  of  simulation 
technique.   The  generation  of  random  numbers  plays  an  important 
role  in  this  method.   First,  the  way  to  generate  random  numbers 
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with  uniform  distribution  is  discussed.   Then  by  a  simple  com- 
putation, they  are  changed  to  random  numbers  with  a  normal  dis- 
tribution.  This  set  of  numbers  will  be  used  to  solve  the  problem 
of  finding  the  probability  density  function  of  a  function  of  the 
normal  random  variables. 

There  are  several  ways  to  generate  random  numbers  with  a 

f>   IP 
uniform  distribution.  '    Most  of  the  available  schemes  use 

the  multiplicative  congruential  method.   This  method  is  con- 
cerned with  generating  sequences  of  nonnegative  integers  by  means 
of  a  congruence  relation,  then  they  are  divided  by  their  mode 
to  get  the  numbers  between  0  and  1. 
The  congruence  relation  is 

Un+1  =  a  Un   (mod  pb)  (3.2.1) 

where  p  denotes  the  number  of  numerals  in  the  number  system 
utilized  by  the  computer  and  b  denotes  the  number  of  digits  in 
a  word.   a  is  a  constant  multiplier  of  the  form  a  =  8t  ±   3 
(t  is  a  positive  integer),  and  Uq,  the  starting  value,  is  an 
odd  integer  to  assure  maximal  periods  for  the  sequences  gener- 
ated by  this  method. 

The  principal  value  of  the  uniform  distribution  for  simu- 
lation techniques  lies  in  its  simplicity  and  in  the  fact  that 
it  can  be  used  to  simulate  random  variables  from  almost  any  kind 
of  probability  distribution.     Therefore  when  the  random  vari- 
able of  uniform  distribution  has  been  generated,  the  random 
variable  of  the  desired  distribution  is  obtained  by  means  of  a 
simple  computation. 
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In  general,  there  are  three  ways  to  generate  normal  random 
numbers  from  the  uniform  ones.   These  are  the  central  limit 
approach,  the  direct  approach,  and  the  fast  procedure.     The 
direct  approach  is  discussed  here  due  to  its  faster  calculation 
and  an  exact  result  compared  with  the  others. 

Given  two  uniform  random  variables  R-,  and  Rp  that  are  inde- 
pendent and  defined  on  the  (0,  1)  interval,  then 

Z1  =    (-2  loge  R-^1/2  cos  2irR2  (3.2.2) 

Z2  =  (-2  loge  R1)1/2  sin  2tiR2  (3.2.3) 

are  two  independent  normal  variables. 

It  is  also  very  easy  to  change  Z   and  Z?  into  a  pair  of 


correlated  normal  variables  with  mean  vector  p,   and  covariance 
matrix  V.   There  exists  a  unique  lower  triangular  matrix  C 


such  that 


X  =  c  z 


where  Z  is  a  standard  uncorrelated  vector  of  random  variables 

as  generated  by  equation  (3.2.2)  and  (3.2.3).   X  is  the  desired 

~>   -*■   -»  m 

correlated  normal  vector,  and  V  =  C  •  C   .   (T  means  transpose.) 

When  two  random  normal  variables  X-,  and  X?  are  given,  their 
means,  variance  and  correlation  coefficient  are  then  defined  as 
EXl  (xx)  =  y.lt    EX2(x2)  =  p.2,  VarXi(x-L)  =  (T^,  Corx^^x^) 

Consequently, 
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V     = 


en, 


r^cr: 


fair. 


c    = 


en, 


P<r, 


cr22 


,VT^ 


(3.2.14.) 


By   equation    (3.2.1j_),    then 


+   H 


*    Txil   -»  lzi 

X  =  =  c 

x2  z2 


<rlZl 


cr2(/Pz1  +  Yi-  ^2  z2) 


M-2 


(3.2.5) 


Equation  (3.2.5)  is  the  equation  to  find  two  correlated  normal 
random  variables.   The  test  of  the  accuracy  of  these  random 
variables  will  be  discussed  in  the  next  section. 

3.3   The  Accuracy  of  the  Monte  Carlo  Method 

The  accuracy  of  the  Monte  Carlo  method  mainly  depends  on 
the  sample  size  of  the  generating  numbers  and  the  validity  of 
randomness.   By  the  law  of  large  numbers,  ^  it  is  shown  that 


p{|x.  -  Pj|<e) 


P1q1 

>  1 Li 


y 


n£ 


-    Xjl  +  Xj2  +  • • •  +  Xjn 

*j 

n 


(3.3.1) 


where   6  is  any  infinitesimal  positive  number, 

p.  is  the  actual  probability  of  occurrence  in  the  1 

J 

interval, 

and    q.  =  1  -  p.  . 
J        j 

n  is  the  sample  size,  and 


th 
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X-jL  =  1  when  the  event  occurs  in  the  jth  interval, 

othe  rwise 
X..  =  0. 

If  the  k  is  the  total  number  of  occurrences  out  of  n  trials 
in  that  interval,  then  X.  =  k/n.   The  value  is  obtained  by  the 

J 

Monte  Carlo  method.   Thus  as  n  — >  oo,  equation  (3.3.1)  yields 
pflXj  -  Pj  |  $  fej— >i  (3-3.2) 

The  above  equation  illustrates  the  probability  of  the  dif- 
ference between  the  value  by  the  Monte  Carlo  method  and  the 
actual  value  which  is  less  than  a  small  quantity.   As  n  approaches 
infinity,  this  probability  becomes  one.   This  means  the  differ- 
ence is  less  than  a  small  quantity  £  when  n  is  sufficiently 
large.   No  matter  how  small  the  £  is,  the  equation  (3.3.2) 
holds.   The  larger  n  the  smaller  6  can  be  expected. 

However,  it  is  impossible  to  increase  the  sample  size  n  to 
infinity.  A   compromise  is  made  considering  the  economic  con- 
sideration of  computer  running  time,  the  size  of  memory  of  com- 
puter and  how  much  accuracy  is  needed. 

Next  question  is  the  validity  of  randomness.   The  statis- 
tical properties  of  random  numbers  generated  by  the  methods  out- 
lined in  the  previous  section  should  coincide  with  the  statis- 
tical properties  of  numbers  generated  by  an  idealized  chance 
device  that  selects  numbers  from  a  certain  interval  independently 
where  each  number  has  a  certain  probability.   Clearly  the  random 
numbers  produced  by  computer  programs  are  not  random  in  this 
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sense,  since  they  are  completely  determined  by  the  starting  'data 
and  have  limited  precision.   But  so  long  as  the  generated  random 
numbers  can  pass  the  set  of  statistical  tests  implied  by  the 
aforementioned  idealized  chance  device,  these  numbers  can  be 

treated  as  truly  random.   There  are  several  tests  used  to  test 

3  cj  TO 
the  validity  of  randomness.  '^,     Here  only  chi-square  test  for 

the  goodness  of  fit  is  discussed. 

With  the  assumption  that  in  each  of  n  independent  trials 

precisely  one  of  r  events,  A-^,    A2,  .  .  . ,  A  must  happen, 

vl»  v2>  *  '  *  >    vn  are  ^e  nurn^ers  °f  successes  out  of  n  in  each 

event  by  the  Monte  Carlo  method.   With  the  actual  probabilities 

r 
plo'  p2o'  *  *  "'  pro  are  numtiers  with  £_     p.   =  1.   In  order  to 

3=1  J 

test  the  hypothesis  p±  =   plQ,  p2  =  p2o,  .  .  .,  pp  =  ppo,  the 
following  statistic  is  considered: 


D2  =  f   (VJ  -  nPJ°)2  =  £     n(pj  -  Pjo> 
3=1     nPj0       j=l      p.o 


(3.3.3) 


It  can  be  shown  that  D2  has  asymptotically  a  chi-square 
distribution  with  r  =  1  degree  of  freedom. 

The  test  of  the  hypothesis  at  the  100  •><  per  cent  signifi- 
cance level  is  obtained  by  choosing  a  number  b  such  that 
P  {"X  >  b|  =  *.   Where  X   has  chi-square  distribution  with  r  -  1 
degrees  of  freedom  and  rejecting  the  hypothesis  if  a  value  of 
Dc  greater  than  b  is  actually  observed.   The  chi-square  value 
for  a  different  degree  of  freedom  can  be  found  in  a  good  mathe- 
matical table. 
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3.1]-     Curve  Fitting  Technique 

Two  kinds  of  fitting  curves  are  used  in  this  paper.   The 
first  kind  corresponds  to  a  polynomial  multiplied  by  an  expo- 
nential term.   This  curve  can  fit  smooth  data  with  a  few  terms 
and  match  its  peak  value  by  the  exponential  term  very  easily. 
The  second  kind  is  a  finite  Fourier  series  of  sine  terms  only. 
It  can  fit  rough  data  faster  than  the  first  kind.   The  fitting 
technique  is  based  on  minimizing  the  least  square  error  between 
the  assumed  curve  and  fitting  data.   The  iterative  method  is 
used  to  find  the  coefficients  of  the  fitting  curve,  with  the 
assumption  that  the  first  kind  of  fitting  curve  is 


f(x)  =  £     arxr  e"Px2  (3.4.1) 

r=o 


Since  the  integration  of  the  density  function  over  its  entire 
range  of  the  variable  is  unity,  the  above  equation  is  modified 
in  order  to  satisfy  this  constraint.   Introducing  a  Lagrangian 
multiplier  A,  the  modified  equation  is 

f(x)  =  f(x)  +  X(   I   f(x)dx  -  1)  (3.1]-. 2) 

'  -  CO 

with  the  assumption  that  g(x^)  is  the  ith  given  datum  of  the 
fitting  set,  and  f(xi)  is  the  ith  datum  of  the  fitted  set.   Then 
the  error  vector  is  E  =  (eA    =  fg(x±)    -  f(xi)|  . 


The  squared  error  is 
where  i  =  1,  2,  .  .  .  n. 


21 


To  minimize  the  squared  error,  a  set  of  equations  are  ob- 


*H 


tained  by  setting  =  0  for  a  particular  a„   and  ~ — ■ —  =  0 


d  e 


2 


r. 


^rn  °      d    \ 


o 


Based  on  these  two  derivatives,  a  matrix  equation  is  written  as 

W*  f •^i}  ■  W   +*  f*"i>]  (3.4.4) 

where  the  notation    [  denotes  a  square  or  rectangular  matrix 
and  ■[  J  denotes  a  column  vector.   The  details  of  the  above 
equation  are  presented  in  Appendix  B. 

For  a  given  p  and  setting  A  =  o,  the  coefficient  vector 
{a^  is  obtained  with  a  value  [a]_1{c}.   Then  by  assuming 
a  X  $   o,  and  the  previous  coefficient  vector  (a J  substituted 
into  |Y(ai)j(  ,  a  new  coefficient  and  vector  can  be  obtained. 

Every  time  a  new  coefficient  and  X  vector  is  obtained  by 
putting  the  old  one  of  a  step  before  into  equation  (3.J+.5)  . 
After  several  steps  it  is  expected  (aA    and  A  will  converge  to 
some  value.   That  is  a  vector  satisfying  the  least  square  error 
fitting  requirement  under  a  preassigned  p.   The  error  ||e||2  can 
be  computed  by  (3.k.3)    immediately.   The  variation  of  this  error 


Note:   When  [a]  is  not  a  square  matrix,  both  sides  of 
equation  of  (3.J±.k)    are  then  multiplied  by  transpose  of  f>1, 
i.e.,  LAJ  .   JAKLAJ  =  B  is  a  square  matrix.   Thus  it  makes  the 
inversion  possible.   So  that  the  least  square  solution  can  be 
obtained   Later  on,  all  the  notation  of  [  ]  is  assumed  as  a 
square  matrix  to  make  the  discussion  of  the  problem  simple 
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function  by  changing  p  can  be  used  as  a  criterion  finding  opti- 
mum p;  i.e.,  whenljEiH   -  \^i-i\   2*   o  and  ||Ei+1  [| 2-  |  [eJI2  >  o. 
there  is  a  local  optimum  Po  such  that  Pi  <  Po  <.   Pi+1,  ||e||  =  f(p) 
in  the  small  neighborhood  about  Po  can  be  found  by  the  tech- 
nique mentioned  above.   The  optimum  Po  is  achieved  by  minimiz- 
ing f(p)  . 

A  Lagrangian  multiplier  X    is  used  to  adjust  the  integrating 
area  of  the  density  function  to  one.   It  can  be  assumed  as  a 
small  value  less  than  0.5,  because  the  result  of  the  density 
functions  by  the  Monte  Carlo  method  does  not  cause  a  large  error. 

Figure  3  is  the  flow  chart  of  the  algorithm  used  in  the 
first  kind  of  curve  fitting.   Appendix  B  has  a  detail  deviation 
of  the  first  kind  curve  fitting  algorithm. 

The  second  kind  of  the  curve  fitting  is  of  the  form 

n 

f(x)  =  £T  ak  sln  kx  (3.I4..6) 

k=o 

Given  f(j  A   x)  is  the  jth  fitting  data,  the  x  is  the  incre- 
ment of  the  fitting  interval.   The  equation  (3.1j..6)  is  written 
as  a  matrix  form 

(f(jAx)j  =  [sin  (jkAx)J  fakJ  (3.^.7) 


with 

2 
n+1 


[sin(jk^  x)]-1  =  ( )  [sin(jkAx)]  , 

therefore 


{akj  =  (  — )  [sin(jkAx)]  (f(jAx))  O.4.8) 
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Pig.  3.   Plow  chart  of  the  first  kind  curve  fitting 
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Fig.  h,.      Flow  chart  of  the  second  kind  of  fitting  curve 


n 


f(x)  =  27  a^  sin  kx. 
k=l 
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By  considering  the  above  equation,  the  inner  products  of 
the  columns  of  the  sine  matrix  with  respect  to  the  data  vector 
f(jAx)   denote  how  much  is  the  contribution  of  the  correspond- 
ing column  to  the  fitting.   A  simple  algorithm  called  the  con- 
strained least  mean  square  solution  is  therefore  suggested.   It 
is  a  method  to  arrange  the  columns  of  the  sine  matrix  according 
to  their  inner  products.   Then,  with  a  given  error  bound,  it  can 
also  adjust  the  mean  squared  error  of  the  fitting  to  within 
this  bound.   Since  this  method  is  designed  to  fit  the  data 
occupied  between  the  interval  0  to  %   and  all  harmonic  terms  of 
sine  function  at  0  and  %   are  zero,  therefore  if  both  or  either 
terminals  of  the  fitting  data  are  nonzero,  it  is  necessary  to 
set  them  to  zero  by  adjusting  the  length  of  the  interval  due  to 
the  convergent  consideration.   Also  for  any  given  interval  T, 
the  scale  factor  T/tc  is  used  to  normalize  this  interval  to  the 
designed  one.   Equation  (3.^.8)  is  thus  rewritten  as 

n         kit 
f(t)  =  H  ak  sin  —  t  (3.IJ..9) 

k=l        T 

Figure  q.  is  the  flow  chart  of  the  algorithm  used  in  the  second 
kind  of  curve  fitting. 

The  methods  discussed  in  this  chapter  will  be  used  to  find 

density  function  of  cos  6-^  in  the  next  chapter. 


CHAPTER  IV 

THE  RESULTS  OP  DENSITY  FUNCTIONS 
BY  THE  MONTE  CARLO  METHOD 

Jj.,1   Introduction 

This  chapter  is  going  to  use  the  methods  of  Chapter  III 
to  find  the  desired  probability  function. 

The  subroutine  RANDU  of  IBM  3&0  is  used  to  generate  two 
sets  of  independent  uniform  random  variables.   These  uniform 
random  variables  are  changed  to  dependent  normal  pairs.   Based 
on  these  pairs,  a  new  random  variable,  that  is,  a  function  of 
normal  random  pairs,  is  obtained.   The  chi-square  test  is  used 
to  check  the  randomness  of  uniform  and  normal  random  variables. 
A  simple  sorting  program  is  run  on  these  random  numbers  to  find 
the  frequency  ratios  of  the  distribution.   After  having  the  fre- 
quency ratio,  the  curve  fitting  technique  is  used  to  find  an 
explicit  expression  for  the  desired  density  function.   In  order 
to  decrease  the  round-off  errors  due  to  the  computation,  all 
of  the  computer  programs  of  this  chapter  use  the  double  preci- 
sion.  In  section  I4..3,  a  comparison  is  made  to  the  densities  of 
S  and  V  by  using  the  Monte  Carlo  method  and  the  analytical  method 
for  the  0~  -   1  and  f  =  O.if,  0.9  cases.   The  densities  of  cos  6-, 
are  obtained  with  explicit  expressions  for  certain  parameters 
in  the  last  section. 
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ij.,2  The  Result  of  the  Generated  Random  Numbers 

First  of  all,  the  IBM  scientific  subroutine  RANDU  is  used 
here  to  generate  a  pair  of  uniform  random  variables.   By  equa- 
tions (3.2.2),  (3.2.3),  and  (3-2.5),  a  transformation  is  made 
on  them  to  pairs  of  dependent  variables. 

This  subroutine  RANDU  is  specific  to  the  system  360  and 

29 
will  produce  2  7  uniform  random  numbers  before  repeating.   The 

seed  a  in  equation  (3-2.1)  is  suggested  to  be  65,539.   21if,358,88l 

and  776,179,721  are  the  two  odd  numbers  used  to  generate  a  pair 

of  uniform  random  variables.   The  sample  size  is  chosen  to  be 

lCn-  throughout  this  paper. 

After  the  uniform  random  variables  are  obtained,  they  are 
sorted  according  to  the  magnitude  of  their  values  and  put  into 
fifty  equal  distance  cells  from  zero  to  one.   Then  the  frequency 
ratio  of  each  cell  is  the  number  of  random  numbers  on  that  cell 
over  the  length  of  cell  multiplied  by  the  total  sample  size. 
Since  the  actual  frequency  ratio  is  unity  everywhere  for  a 
uniform  distribution,  chi-square  values  of  these  pairs  of  uni- 
form random  numbers  are  calculated  by  the  definition  of  section 
3-3.   They  are  6£.lj.9  and  l±3.l\.6   respectively.   For  the  case  of 
k9   degrees  of  freedom,  P  \%2   >  77-7331+]  =  0.005  is  the  value. 
Because  the  above  two  values  are  less  than  77.7331+,  these  two 
sets  of  uniform  random  numbers  are  not  rejected  for  0.5  per 
cent  significance  level. 

By  equations  (3-22)  and  (3-23)  two  independent  normal 
variables  with  zero  means  and  standard  deviations  of  one  are 
obtained.   Although  the  range  of  these  random  numbers  is  between 
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positive  and  negative  infinities,  it  is  cut  down  to  a  value  from 
+3.5  to  -3.5  (i.e.,  =  3.5  <1T~).   That  is  reasonable  since  the 
probability  that  the  normal  random  numbers  are  outside  the  afore- 
mentioned range  is  less  than  O.OI4.7.   For  finding  the  chi-square 
of  the  generated  normal  random  variables,  the  same  procedure  is 
done  as  uniform  random  variable  case.   The  only  difference  is 
its  range  of  the  variable.   In  this  case,  the  chi-square  values 
are  7&.20  and  61|. 20,  respectively.   They  are  also  not  rejected 
for  the  0.5  per  cent  significance  level.   Figure  5  shows  a. histo- 
gram of  one  of  the  generated  normal  density  function.   The  dot 
points  express  the  value  in  the  exact  normal  density  function. 
By  using  a  curve  fitting  method  to  fit  the  histogram,  those 
irregular  deviations  caused  by  errors  are  smoothed  out. 

lj..3   Comparisons  of  the  Results  by  the  Monte  Carlo  Method 

and  the  Analytical  Method 

The  density  functions  of  S  and  V  are  obtained  in  this  sec- 
tion by  the  Monte  Carlo  method.   They  are  compared  with  the 
results  obtained  by  analytical  methods  in  section  2.3. 

The  two  sets  of  independent  normal  variables  which  are 
found  in  the  last  section  are  changed  to  dependent  ones  with  a 
given  standard  deviation  vector  and  correlation  coefficient. 
The  mean  vector  is  assumed  to  be  zero  and  the  two  random  vari- 
ables X  and  Y  have  the  same  standard  deviation  throughout  this 
paper.   The  random  variables,  S  and  V,  are  calculated  by  their 
functional  relationship  with  the  generated  random  variables. 
The  range  of  V  is  from  zero  to  one  and  S  is  cut  down  from  one 
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Fig.  5«   Histogram  of  the  generated  normal  distribution 

of  zero  mean  and  standard  deviation  one 

by  the  Monte  Carlo  method. 


30 

to  5-0.   By  the  same  method  as  used  in  the  last  section,  they 
are  sorted  according  to  the  magnitude  of  their  values  into  fifty 
equal  distance  cells  in  the  range  of  the  variable.   The  frequency 
ratios  of  each  cell  is  calculated  thereafter.   In  other  words, 
the  histograms  of  S  and  V  are  found.   The  first  kind  of  fitting 
curve  is  then  used  to  fit  these  histograms.   The  coefficients 
of  the  fitted  equations  S  and  V  for  Q~  -   1  and  P  -   O.lj.,  0.9 
are  listed  in  Table  1.   Actually,  the  fitted  equations  express 
the  density  functions  of  V  and  S  with  the  given  parameters. 
Their  curves  are  plotted  in  Pig.  6  to  Pig.  9.   They  are  compared 
with  the  corresponding  curves  obtained  by  the  Monte  Carlo 
method.   These  comparisons  show  that  fairly  close  results  are 
obtained  by  both  methods.   Because  the  inverse  quantity  might 
cause  larger  round-off  error,  the  densities  of  V  produce  more 
error  than  the  densities  of  S  by  the  Monte  Carlo  method.   How- 
ever, the  errors  are  not  over  +3  per  cent.   These  results  verify 
that  the  Monte  Carlo  method  is  a  good  valid  approximation.   By 
using  the  Monte  Carlo  method  the  probability  density  of  cos  0-, 
will  be  obtained  in  the  next  section. 

W     The  Density  Function  of  Cos  Q-^   by  the 

Monte  Carlo  Method 

Although  by  using  analytical  methods  to  find  an  algebraic 
solution  of  the  density  of  the  random  variable  cos  6-i  is  very 
hard,  the  Monte  Carlo  method  can  do  this  job  numerically  without 
too  many  difficulties. 
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Table  1.   Coefficients  of  fitting  curves  of 
S  and.  V  with  QT  =  1 


f(x) 


-I 


r=0 


arxre-P(x-Pl} 


f=   0.4 


P=   0.9 


P=  o.k 


0.9 


Pi 

0.705 

0.0 

0.0 

0.0 

p 

11.7781; 

1.5137 

0.779 

0.552 

a0 

-18.00380 

1.37614 

-7.37264 

42.52669 

al 

132. 42396 

-22.27794 

17.24265 

-86.40062 

a2 

-292.89944 

104.01710 

-9.15744 

69.11961 

a3 

272.25791 

-157.98395 

0.50329 

-24.80488 

a4 

-86.77027 

85.05811 

1.01621 

3.45576 
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0.0   . 


By  the  Monte  Carlo  method 


By  the  analytical  method 


Pig.  7.   Density  functions  of  S  with  (7"  = 

f>=   0.9. 
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By  the  same  method  as  section  I+..2,  the  new  random  variable 
cos  0-,  is  obtained  by  its  functional  relationship  to  X  and  Y  in 
the  cases  CT   =  0.1,  3,  100,  f>=   0.1].,  0.9,  and  0  =  30°,  k$° ,  60°. 
Frequency  ratios  of  each  cos  0-,  with  different  parameters  are 
calculated.   The  first  kind  of  curve  fitting  technique  is  then 
used  for  Q~   =  0.1  case.   However,  for  CT  >  1,  the  second  kind 
of  curve  fitting  is  used  due  to  the  histograms  changing  so 
rapidly.   In  Table  2  to  Table  7  are  listed  all  the  coefficients 
of  the  fitted  curves.   Pairs  of  density  curves  of  cos  0-^  with 
equal  o~   and  0  but  different  P   are  plotted  from  Pig.  10  to 
Fig.  18.   Since  the  total  areas  under  the  densities  are  unity, 
the  areas  in  Fig.  l£  are  calculated  by  the  trapezoidal  method 
as  a  valid  check.   They  are  0.951  for  P  =   0.9,  and  1.0098  for 
f  =  O.ij.,  respectively. 

From  these  curves  and  tables,  observations  can  be  made  with 
respect  to  changing  of  parameters  U~,     P  ,  and  0.   The  range  of 
the  random  variable  cos  ©^  is  between  +1  and  -1.   For  CT   =  0.1 
and  the  different  r    cases,  the  curve  occupies  only  the  posi- 
tive axis  with  one  peak  value  (Figs.  10-12)  .   When  C7~  increases, 
it  extends  to  the  negative  axis  with  two  major  peak  values  in 
the  positive  and  negative  sides  respectively  (Figs.  13-18). 
The  amplitudes  of  these  densities  are  also  affected  by  CT .      If 
the  standard  deviation  is  smaller,  then  the  amplitude  of  the 
peak  is  larger. 

The  correlation  coefficient  does  not  cause  a  significant 
difference  of  the  densities  for  the  CT  =   0.1  case  (Figs.  10-12). 
It  does  not  affect  the  density  greatly  when  CT   increases.   The 
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Fig.  10.   Density  functions  of  cos  9]_  with 

o-  =  o.i,  e  =  30°. 
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f(cos  6-j^) 


0.440 


0.706 


Pig.  11.   Density  functions  of  cos  0-. 
with  <T  =  0.1,  0  =  45°. 
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f(cos  e1) 


0.116 


Fig.  12.   Density  functions  of  cos  0-|_ 
with  (T  =   0.1,  9  =  60°. 
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larger  P  produces  a  higher  peak.   The  range  of  cos  6^  also 
Increases  when  P   decreases  (Pigs.  13-18) . 

When  cr  is  equal  to  0.1,  the  angle  9  exactly  locates  the 
point  where  the  peak  value  of  the  densities  occurs.   That  is 
the  value  of  cos  0  (Figs.  10-12).   However,  when  <J~   is  equal  to 
3  and  100,  it  is  hard  to  see  where  the  peak  value  does  occur. 
A  conservative  estimate  is  around  the  value  of  cos  0  for  the 
major  peak  of  the  positive  axis.   The  major  peak  of  the  negative 
axis  is  closer  to  the  origin  than  the  peak  of  the  positive  axis. 
This  estimate  is  better  for  the  cr  =  3  case  than  the  0~  =  100 
case  (Pigs.  13-18) . 

As  a  physical  interpretation,  when  a~    is  large  sharp  slopes 
in  the  reflecting  body  are  expected.   As  example  <T~   =  100,  over 
90  per  cent  of  the  surface  slopes  are  between  1+5°  and  87  • 
Especially  when  0  is  very  large,  the  most  values  of  cos  0-j_  are 
either  in  the  same  direction  as  a^  or  in  the  opposite  direction 
(i.e.,  cos  0  — >  ±1) .   When  0  decreases,  a  considerable  number 
of  cos  0-,  values  approach  0-,  =  90°  (i.e.,  cos  0-^  =  0).   These 
phenomena  are  clearly  shown  in  Fig.  19.   They  are  actually 
observed  in  Figs.  13  to  18. 

From  the  tables,  when  0~    =0.1,  the  densities  of  cos  0-^ 
are  approximately  normal  distributions  with  mean  0.5  for  0  =  60° 
as  shown  (Table  3).   The  densities  are  centered  at  0.707  and 
0.866  for  0  =  \\S°   and  0  =  60°,  respectively,  with  the  same 
standard  deviation  0.1  (Table  2,  Table  \\.)  .   They  are  somewhat 
like  Rice  distribution  for  the  latter  two  cases.   The  existence 
of  higher  order  harmonic  terms  in  the  fitting  curve  of  cos  0n 
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for  the  cases  CT"  =  3  and  CT  =  100  means  more  roughness  in  the 
data  as  given  by  the  Monte  Carlo  method.   Theoretically,  such 
roughness  cannot  occur.   This  is  an  accumulated  error  due  to 
the  curve  fitting  and  the  computer  simulation.   The  error  of 
curve  fitting  is  not  so  great  because  it  uses  the  minimum  least 
square  sense.   The  computer  simulation  error  is  actually  the 
origin  of  the  error.   By  increasing  the  sample  size,  a  better 
result  is  expected.   However,  the  computer  program  used  in  this 
paper  has  to  use  a  double  precision  and  pairs  of  generated 
random  numbers  have  to  be  stored.   Therefore  by  using  IBM 
360/50  with  the  double  precision  and  the  sample  size  10,000 
(for  pairs,  i.e.,  20,000),  the  size  cannot  be  increased  sig- 
nificantly.  But  to  avoid  the  problem  of  the  limited  memory  of 
the  computer,  one  way  is  to  store  the  generated  random  numbers 
on  tapes.   This  is  a  recommended  suggestion  for  the  future 
research. 


Table  2.   Coefficients  of  fitting  curves  of  cos  9-p 
with  ir  s  o.l,  9  =  30°, 

f(x)  =  £  arx  e  HV   K1/  . 
r=0 


k9 


<r  =   o.l,  e  =  30° 


/*  =  0.I4.  P  =  0.9 


Pi 


0.866  0.866 

p                  192.247  2014..  370 

a0                3363.18032  7859.27189 

a1                                  -II844.72088  -27538.79014 

a2              13917. 81+435  32136.8961 

a,               -5443.13619  -12477.08512 


Table  3.   Coefficients  of  fitting  curves  of  cos  8^ 
with  <T  =   o.l,  6  =  60,  f(x)  =  Ae"p(x_pl) 


cr  =  o.l,  6  =  60 


P  =   0.4  p  =   0.9 


Px  0.5  0.5 

p  69.928  67.889 

A  4.72168  4.65298 
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Table  1]..   Coefficients  of  fitting  curves  of  cos  6-1 
with  <T-   0.1,  6  =  \\$° , 

*t    \         rr-     r  -p(x-P"i  ) 
f(x)  =  £__  arx  e  ^v   ^J-'  . 

r=0 


cr  =  o.i,  e  =  k$° 


P  =  0.24-  /*  =  0.9 


Pi  0.707  0.707 

p  101.859  120.00 

aQ  265.329  692.26728 

a1  -1130.96678  -2877.5886i|. 

a2  1635.07361  3998.563^2 

83  -78if.38536  -I8I4I.2388I4 
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Table  5>.   Coefficients  of  fitting  curves  of  cos  6-, 
with  <T  =   3,  6  =  30°. 


(°=o.k 


P=   0.9 


Order 
No. 


Coefficient 


Order 
No. 


Coefficient 


1 

0.99702 

3 

0.29310 

k 

0.10095 

6 

0.0701+9 

2 

-O.069I4.O 

10 

-0.06015 

12 

-0.03525 

18 

-0.03381+ 

27 

0.033i|-9 

5 

-0.03190 

13 

-0.021+20 

1 

0.98888 

3 

0.1+1367 

k 

>     0.36381^ 

7 

0.22l].71 

2 

-0.15871 

Terms 
of 
fitting 


11 


Terms 
of 
fitting 
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Table  6. 


Order 
No. 


Coefficients  of  fitting  curves  of  cos  9-. 
with  or   =  3,  e  =  ij.50.  x 


e=  o.k 


Coefficient 


l°=   0.9 


Order 
No. 


Coefficient 


1 

0.77251^ 

1 

0.76162 

3 

0.39098 

3 

0.60205 

5 

0.18252 

6 

O.Uj.228 

7 

O.08083 

2 

-0.l30i|.7 

k 

-0.079098 

9 

0.10959 

2 

-0. 07^57 

10 

-0.10710 

6 

-O.Olj.739 

k 

0.05369 

21 

O.02976 

5 

0.08159 

1$ 

-0.03256 

8 

-0.05959 

11 

-0.03134 

12 

-0. 014.341 

20 

0.03166 

20 

0.03795 

38 

0.08162 

33 

-0.02759 

28 

-0.02775 

Terms 
of 
fitting 


111- 


Terms 
of 
fitting 


11 


Table  7-   Coefficient  of  fitting  curve  of  cos  9^ 
wither  =  3,  0  =  60°. 
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Order 
No. 


P=  o.k 


Coefficient 


/"  =  0.9 


Order 
No. 


Coefficient 


1 

0.61517 

3 

0.34107 

5 

0 .  20143 

7 

0.16137 

9 

0.13458 

16 

-0.08781 

12 

-0.08355 

11 

0.07622 

10 

0.07668 

13 

0.07392 

8 

-0.07460 

18 

-0.07502 

6 

-0.06877 

14 

-0.06860 

15 

0.06057 

17 

0.04955 

22 

-0.05024 

Range 


Terms 
of 
fitting 


-0.8 


1.0 


17 


1 
3 
5 

14 
12 
22 

2 
16 


Range 


Terms 
of 
fitting 


0.62455 

0.60531 

0.26235 

-0.05642 

-0.04905 
-0.04840 
-O.03867 
-0.03540 


■0.7 


1.0 


8 
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Table  8.   Coefficients  of  fitting  curves  of  cos  9-, 

with  cr  =  ioo,  e  =  30°. 


e- 

=  0.4 

P- 

0.9 

Order: 
No.  : 

Coef-    : 
ficient  : 

Order: 

No.  : 

Coef- 
ficient 

Order: 

:  No.  : 

Coef- 
ficient 

:  Order: 

:   No.  : 

Coef- 
ficient 

1 

1.0091*4 

11+ 

0 . 14006 

3 

1.11162 

8 

0.05313 

3 

0.68269 

26 

0.09201 

1 

1.01148 

22 

0.05175 

5 

O.48581 

12 

0.11103 

5 

0.41252 

24 

0.05059 

7 

0.37063 

15 

0.10281 

6 

0.19378 

1+1 

0.04973 

9 

0.29316 

28 

0.05442 

1+ 

0.13130 

11 

0.21800 

10 

0.08699 

13 

0.09265 

18 

0.20250 

31 

0.11149 

11 

0.07872 

16 

o.i9i5ov 

3^ 

0.11179 

16 

0.07023 

20 

0.18180 

32 

0.02277 

15 

0.06885 

f 

22 

0.16734 

k2> 

0.09621 

18 

0.06152 

24 

0.13745 

10 

-0.05963 

IB 

0. 11+415 

11+ 

0.05315 

Rang 

e    :    -0 

.48  — 

0.58     r 

Ran 

ge    : 

-0.48  — 

0.58 

Terms 
of 
fitting 

22 

Terms 
of 
fitting 

16 
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Table  9.   Coefficients  of  fitting  curves  of  cos  Q^_ 
with  <T    =   100,  0  =  k$°. 


/° 

=  0.1| 

P 

=  0.9 

Order: 
No.    : 

Coef-         : 
ficient    : 

Order: 
No.     : 

Coef- 
ficient 

Order: 

No.     : 

Coef- 
ficient 

:    Order: 

:      No.    : 

Coef- 
ficient 

1 

0.73180 

kk 

-O.I3236 

3 

O.80616 

27 

-0.061j.88 

3 

O.ij.9717 

16 

0.13527 

1 

0.71i|68 

k2 

-0.05182 

5 

0.37588 

k2 

-0.11895 

5 

0.32176 

ko 

-0.05171 

7 

0.2902l| 

2k 

0.12506 

6 

0.09813 

16 

0.01^123 

9 

0.23121 

36 

-0.12ij45 

13 

0.07023 

25 

-0.05101 

11 

0.16809 

Ik 

0.12319 

kh 

-0.07021 

2i+ 

0.03620 

22 

0.15313 

27 

-0.09990 

k 

O.06139 

18 

0 . 15035 

15 

0.05217 

11 

0 .  05lif7 

38 

-0.1i|206 

29 

-O.09i4.i3 

21 

-O.06185 

13 

0.09973 

26 

0.11080 

29 

-0.05998 

20 

0.1^376 

23 

-0.05352 

U-0 

-0.12922 

ks 

-0.05910 

Rang 

e        : 

0.68  — 

O.78           : 

Ran 

ge        : 

-0.68  — 

O.78 

Terms 
of 
fitting 

22 

Terms 
of 
fitting 

i8 
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Table  10.   Coefficients  of  fitting  curves  of  cos  0-j_ 
with  cr  =  100,  0  =  60°. 


r- 

:  O.I4 

• 
« 

• 

P- 

0.9 

Order: 
No.  : 

Coef-    : 
ficient  : 

Order: 

No.  : 

Coef-    : 
ficient  : 

Order: 
No.  : 

Coef- 
ficient 

:  Order: 
:   No.  : 

Coef- 
ficient 

1 

0.57913 

27 

-0.07782 

3 

0.61967 

36 

-O.O39I46 

3 

0.38075 

38 

-0.07357 

1 

0.56811 

I4O 

-O.OI4306 

5 

0.29065 

26 

0.07253 

5 

0.2l(.5l7 

7 

0.22918 

16 

0.06938 

6 

0.09974 

9 

O.I6796 

¥> 

-0.06857 

13 

O.0253I4. 

11 

0.12681 

29 

-0.06596 

25 

-O.069I4.3 

22 

0.10099 

36 

-0.05705 

23 

-0.07201 

2k 

0.09531J- 

ki 

0.0560i| 

k 

0.05738 

20 

0.09171 

1+8 

-0.05523 

21 

-0.07097 

13 

O.090I4.5 

ki 

0.05337 

kh 

-O.OI4.892 

kh 

-0.08523 

ik 

0.05113 

11 

0.01302 

18 

0.08226 

28 

0.024.7814. 

38 

-O.Ol4J4.Ol 

ko 

-0.08090 

15 

O.Olj.715 

2k 

0.05995 

2$ 

-O.O80I4J4. 

21 

-O.Ol4.657 

18 

0.06820 

23 

-0.07956 

31 

-O.OI1.382 

k2 

-  0.01+563 

k2 

-0.07765 

kl 

0 .  OI4.3I4.8 

10 

-0.01168 

Rang 

e    :     -0 

.86  — 

0.98    : 

Range    : 

-0.86  — 

0.98 

Terras 
of 
fitting 

32 

Terras 
of 
fitting 

18 

CHAPTER  V 

CONCLUSIONS  AND  RECOMMENDATIONS 

5.1  Conclusions 

The  goal  of  finding  the  probability  density  of  cos  Q1   which 
is  a  function  of  two  correlated  random  variables  X  and  Y  does 
not  succeed  by  using  the  analytical  methods.   Therefore  the  Monte 
Carlo  method  is  used.   When  this  method  is  combined  with  the 
curve  fitting  techniques,  approximate  solutions  are  obtained  in 
a  closed  form  for  different  parameters  a~ ,  f> ,    ©.   Three  con- 
clusions are  made  to  these  results. 

1.  When  the  standard  deviations  of  the  two  random  variables 
X  and  Y  are  less  than  unity,  the  densities  of  cos  0^  are 
slightly  affected  by  the  changing  of  the  correlation 
coefficient.   They  are  approximately  Rice  distributions 
centered  at  the  value  cos  9. 

2.  When  the  standard  deviations  of  X  and  Y  are  larger  than 
unity,  the  effect  on  the  densities  due  to  increasing 
the  correlation  coefficient  is  significant.   Two  major 
peak  values  of  the  densities  occur  in  positive  and  nega- 
tive sides,  respectively. 

3.  The  larger  the  standard  deviation  of  X  and  Y,  the  larger 
the  occupied  range  of  the  density  of  cos  0j_  is.   However, 
the  range  is  limited  between  -1  to  1. 
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£.2   Recommendations 

1.  One  way  to  find  the  density  of  cos  61   by  analytical 
method  might  be  very  interesting.   Because  the  density 
of  cos  6-,  can  be  obtained  by  the  densities  of 


Y/Vl  +  X2  +  Y2  which  is  the  derivative  of  V 1  +  X2  +  Y2 
with  respect  to  Y.   In  other  words,  with  given  prob- 
ability densities  of  f(X,  Y) ,  the  study  of  finding  the 

a  d 

densities  of  —  f (X,  Y)  or  —  f (X,  Y)  is  suggested. 

Hopefully,  this  approach  will  not  be  tied  up  by  hard 
integrations . 

2.  In  order  to  solve  the  problem  of  the  limited  memory  of 
the  computer,  storing  the  generated  random  numbers  on 
tapes  is  suggested.   Thus  the  sample  size  of  the  Monte 
Carlo  method  can  be  increased  greatly.   A  better  accuracy 
is  then  expected  by  this  method  in  finding  the  densities 
of  cos  0-^. 

3.  With  a  larger  sample  size,. the  use  of  the  Monte  Carlo 
method  to  find  the  probability  density  functions  of  the 
reflection  coefficients  V   (0-^  and  Vnn(01)  is  recom- 
mended, where 


n-L2  cos  6]_  -•Vn22  +  cos2  61 
n-j^2  cos  0-^  +*y  n22  +  cos2  0, 

no  cos  0-l  -  y   n22  +  cos2  Q-^ 
n-,  cos  0,  +Vn  2  +  cos2  eT 

nl>  n2»  n1  are  deterministic  constants. 


vpp(si)  ■ 


vnn(ei)  = 
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ij..  The  densities  of  cos  6-^  change  from  Rice  distribution 
to  two-peak  distributions  when  cr  are  altered  from 
0.1  to  3.   An  interesting  problem  is  to  find  the  value 
of  0~   when  this  changing  exactly  occurs. 
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APPENDIX  A 
A.l  The  Derivation  of  the  Probability  Density  of  cos  6-, 

<a   +  bj>   =   a  +  b<y>  =  a  (A.  1.1) 

<y=r=i>=  <WTT77> +  h<^  77  +  *>     u'1;2) 

By  using  polar   coordinates,    with  c   =   2CT2(l   -  ffi) ,    then 

/         y         v  (     f2%   *Y  1-  f2   r2    sin  6  /     r2  \ 

\T)       p      2/=    /  w  9         exP (l-^sin   2e)}d9dr 

Ni/l+x2+yV        /0    /0  rccYl  +  r2  (      c  7 


f*  YTT^r2  r2     ^  ,r2/*  . 

exp    (-  — )   /  sin  6    exp/  sin   26  f    d6dr 


b      nc-yi+r2  c     h 


where  the  second  integration  of  the  above  equation  is 

r2%  ,         pv2  ,<*>  .  pv2 

sin0jlo( )    +2[Y.      (-I)*  I2k+1    ( )    sin(2(2k+l)©) 

/0         L    c       V  k=0  x        c 


0£> 

+ 


JL     (-Dk  I2k( )  cos  k   ke)[  de  =  °  (A. 1.3) 

k=l  c  f  ) 


The  mean  value  of  cos  Q-,    i; 


/     a  +  by     n  =    y        1         \ 

Xi/l   +    x2    +   y2/         a  \f  1  +   X2    +  y2/ 

f      rYl-/72  r2  r2/> 

=    2a    /  V 5-    exp    ( }    Io    ( )dr  (A. 1.^) 

/0      c-Tl+r2  c  c 
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Assume      z   = 


71  +  x2   +  y2 


F7(z)    = 


>-C*2        -  DO 


2CT2  Yl    -  p2    71 


(x2   +  y2   -    2/^xy) 

• -r, — r^ r  dxdy 

2Cr2(i   _^^) 


By   taking   the    derivative  with  respect   to   z,    it  yields 


Ma)    = 


'-C*£>   C7l(l-Z2)^/2 


/"      (x2+z2)>)            /  2  /^xz  -fl+x2  | 
exp  -  -  (  exp  J  ,  \ 

c(i-z2)J       (  c-yi  -  z    y 


dx 


(A. 1.5) 


For 


^=  0   case,    c   =   2CP 


-D^ 


f7i(z) 


-f  1   +   x2 


(x2+z2) 


exp 


'-oo  2oicr2(l-z2)3/2  2cr2(l-z2) 


dx 


(A. 1.6) 


Equations  (A.l.Ij.),  (A. 1.5),  (A. 1.6)  have  not  been  solved  to 
obtain  the  solutions  which  can  be  expressed  in  algebraic  form 
or  with  well-known  series  expressions.   The  moment  theorem 
states  that 


<g   (jW)n 

<P  (w)  =  2^    ] 


n=0    n 


n 


th  T 

where  mn  is  the  n   moments.   If  they  are  known,  then  (f)  (w) 


can  be  obtained.   Then, 

f(x)  =  —  1    I  (w)  e"Jwx  dx 
2%    /-OoJ' 
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<cosne^> 


-<( 


r   a  +  by 


)>■<■ 


n(n-l) 
an+nban-1y+  b2an_2y2+. . .+bnyn 


21 


(1   +   X2   +   y2)n/2 


n(n-l) 


nban_1y 


:  ^  ( 1+ x2+y2 ) n/2/  +  \ i+x2+y2 ) n/2'      ^ (1  +   x2   +  y2) n/2 


;>< 


2  J 


b2an-2y2 


;> 


n  ,„n 


<- 


b"  y 


(1   +   x2   +  y2)V2 


> 


-> 


The  average  value  of  the  terms  with  odd  power  of  y  vanish  in 
the  above  equation.   The  terms  with  even  powers  of  y,  by  using 
polar  coordinates  and  assuming  n  =  2m,  then 

2™    Vi.^2 


,2m 


<- 


(l+x2+y2) 


a 2oi      (r   sin  9) 
(l+r2)n/2  ^c 


exp 


(1   -  /°  sin   29) 


J  I 


y  d6dr 


where   c   is    defined   as   above,    it  yields 
r2m 


N(l+x2+Y2)n/2/        y0 


r2m    Ap[Tp^ 


2mTtc(l+r2)n/2 


2/       Z9^ 
"r  /C  T   (i— )  dr   (A. 1.7) 


The  above  equation  has  not  been  solved  to  find  a  solution 
which  can  be  expressed  algebraically  or  with  a  well-known 
series  expression. 
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APPENDIX  B 


B.l  The  Derivation  of  the  First  Kind  of  Curve  Fitting  Algorithm 
By  substituting  equation  (3.I4..I)  into  (3.1^.2),  then 


m  o      /     m 


f(x)  =  (£*  arxr)  e"Px  +  \(  )  £  arxr  e"Px   dz  -  I) 

r=0  /-c*>  r=0 

the  square  error  is 

E  |2  =  I!  fg(Xi)  -  ^1  x±r  e"Pxi^  +  \(£  ar[   x^'P*  dx-^J 
i=l  u        r=0  r=0    '-<& 


2 


2 


=  0,  then 


3  \ 

.C*S> 


m  f —       2 

I  ar      xr  e_Px  dx  =  1  (B.l.l) 

r=0    /-« 


<$o 


—  =  0,  with  equation  (B.l),  then 


2% 

i=l     r=0     r   1      x  i=l  x      Loo 


dx 


n         m  2  r         t.  9  n  ? 

=  X    r       ZT     arx.r   e"Pxi  xP0  e-Px'   dx  +    £     x.r0e"Pxi      g(x,) 

i=l     r=0     r  x  Loc  i=l  1 

(B.l. 2) 

where   rQ  =  0,    1,    2,    .     .    . ,    m. 

Equations  (B.l)  and  (B.2)  are  combined  and  written  in  a 
matrix  form. 


i  ffcl. 


;}+x(Y(a.)j 


where 
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(B.1.3) 


W 


a, 


M- 


m  y 

fZTxi0   e"Pxi2   g(xi)' 
^x^   e"Pxi2   g(Xl) 


rXjLm  e"pXi    g(Xi) 


/j-x°   e-P^   dx^ 


fx     e"P^   dx 
Y(ai)    =    {     . 


n  m 


(x™  e"px2   dx 

(  •     ) 


z.   zr  a-x--r  e_pxi' 


i=l      r=0 


r   1 


I 

A21      I       A22 

A  is  a  (m+2)  x  (m+2)  square  matrix 

Ai:l  is  (m+1)  x  (m+1),  A,p  is  (m+1)  x  1,  Ap,  is  1  x  (rat+1), 

App  is  0. 


"See  the  matrix  notation  in  the  footnote  of  section  3.1±. 
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11 


£    x±0x.0e-2pXi2  f         1  Oe-2pXi2     £     x^x.V2?*!2 
i=l  X   X  1=1  X     X  1=1  1  X 


£    x.°x.V2Pxi2,  r  x^x.V^i2,...,^  x.-x.^e 


1  =  1 


i=l 


1=1 

n 

i=l 


m   1  -2px,-' 
x.  « 
l   l 


Jl         0  m  -2PX-;2   "    x    -2pxi2 
2.  x  x  e      ,  £_  x.  x.  e     x  ,  . 


1=1 


l   l 


12 


£    g(x±)  fx° 

i=l    x  ' 

n        ,  1 
T  g(x±)   x 
1=1       i 


1=1 


e"Px  dx 


-  Px   j 
e     dx 


i   i 


f    x.\.me-2pxi2 
1=1  x      x 


n        / 
1=1      / 


xm  e"Px  dx 


l21 


=  [(x°e-Px2  dx,  /x^-P^  dx,  .  .  .,  fA-P*2  dx] 


B.2   The  Algorithm  for  Finding  an  Optimum  P 

With  the  assumption  that  the  exact  error  curve  is 

Ex(p)  =  aQ  +  axp  +  a2p2 

The  calculated  error  value  is  Ed(p.),  then  the  square 
error  between  these  two  values  is 


IMI2  =  f  [ei(pi)  -Ed(pi)]: 
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?  ei 


3  a, 


=  0  ,  where  r  =  0,  1,  2,    then 


£  f2  E-lCpj;)  -  2  Ed(Pl)l  =  0 
i=0  L  ■  J 

f  [2  P^CPi)  -  2  PlEd(Pl)]  =  0 

1=0  ■" 

r  [2  P^E^p.)  -  2  Pi2Ed(Pl)]=0 


The  above  three  equations  are  combined  and  written  in  a 
matrix  form. 


i=l 


n  n 


n 


r  1    r  Pi    i^ 

i=l  i=l  i=l 


n 


n  n 

r  Pi  x:  Pi2  t  Pl3 


i=l  i=l 


i=l 


n  n     n  n 

r  Pi2  r  p^  r  Pi 

1=1  1=1  1=1 


'0 


i=l 


=  V 


,     2 

V     / 


n 


i=l 

n 


1    1 


(B.2.1) 


Vli  p— 


'0 
The  coefficient  vector  •/  a-,  I  can  be  obtained  by  the  equa- 

a2 


tion  (B.2.1).   The  oPtimum  P  is  then  p>  = 


-a- 


2a' 
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APPENDIX  C 
C.l  PROGRAM  FOR  GENERATION  UNIFORM  RANOCM  NUMBERS 


2 

6 

8 

40 

42 

78 


85 
86 


10 


15 

11 


18 


IMPLIC 

DIMENS 

01  MENS 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

READ  ( 

00  3  K 

READ  ( 

DO  3  I 

IY=IX* 

IF  (IV 

IY=IY+ 

YFL=IY 

YFL=YF 

IX=IY 

R(K,I) 

CM=M 

BN=NN 

B(1)=C 

DO  10 

Bin=e 

DO  9  K 
DO  9  I 
C(K, I  ) 
DO  11 
DO  11 
H=0.C2 
DO  15 
IF  (R( 
H=H+G. 
C(K,  I  ) 
DO  18 
'CD(K)= 
AM(K)= 
SM(K)= 
DC  18 
FR(K,  I 
DD(K)= 
AM(K)= 
SM(K)= 
DO  30 
IF  (Af 
DM(K)= 


IT  REAL*8(A-HfG-Z) » INTEGER  I I-N) 
ION  R(2,ICCCG),C(2,5G),B(5C),FR(2,50) 
I CM  AM (2) , DH (2 ) f  SM 1 2  J  t  DD 1 2 ) 
(2110) 

(14X»*0X*t6Xv  'URSeX,  'FRISIX,  »FR2'  ) 
(/12X,F5.2,3X,F7.4,3X,F7,4,3X,r7.<,) 

l/i2X,«MEAN=«  iF9.7,3X,«   VARIAMCE  = 


»,F9.7) 


■  I////12X, 'UNIFORM 
(/12X,'CHI  SCUARE 
1,2)  NN.M 
=  1,2 
1,2)  IX 
=  1,M 
65539 

)  85,86,86 
2147463647+1 

L*.4656613E-9 

=  YFL 


.0 

1=2, NN 
(I-D+0.02 
=  1.2 

=  1,NN 
=  0.0 
K=l,2 
J=l,M 


DISTRiBUl IGN 
VALUE=',F9.4) 


li) 


1  =  1 

K,  J 

C2 

=  C( 

K=l 

CO 

CO 

CO 

1  =  1 

)=C 

cot 

(C( 
(C( 
K=l 
(K) 
SM( 


VNN 

).LE.H)  GO  TO  11 

K,  I  )+1.0 
,2 


•  NN 

(K,I  )*BN/CM 

K)+( (FR(K»I )-l.G)**2)*2QQ 

K,  t)*(B(  I)  +  ,01)/CM)+AfMK) 

K,  I  )*(  {  C<  {  I  »  +  .0i)**2)/CM)+SM(K) 

.2 

.EC.C. )  GO  TO  31 

K)-(AM(K)**2) 
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31 

30 


60 


GO  TO  30 
DM(K)=SM(KJ 

CONTINUE 
U  =  l. 


WRITE 
WRITE 
DO  6C 
WRITE 
WRITE 


(3,6) 

(3,8)  (B(I),U,FR(l,I)tFR{2,I),I=l,NNJ 

1  =  1,2 
(3,42)  I 
(3,40) 


WRITE(3,78) 

STOP 

END 


AM C lit DMCIJ 

CD(  I) 
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C.2  PROGRAM  FOR  GENERATING  INDEPENDENT  NORMAL  RANDOM  NUMBERS 

IMPLICIT  RE AL*8(A-H,0-Z), INTEGER! I-N) 

DIMENSION  Z(2,100CO),C(2,50),B(51),FRi2,50),D(5l) 

DIMENSION  AM (2) , DM ( 2 ) , SM (2 ) , CD ( 2 ) ,G ( 5C ) 

2  FORMAT  (2110) 

6  FORMAT  (14X,«CX»,6X, »NR«, 8X, • FR1 • ,7X, « FR2* ) 

8  FORMAT  l/12X,F5.2, 3X, F7.4, 3X ,F7.4 , 3X, F7. 4 } 

40  FORMAT  (/12X,«MEAN=» ,F9.7,3X,'   VARIANCE      =«,F9.7) 
42  FORMAT  (////12X, 'NORMAL  DISTRIBUTION  »,I1) 
78  FORMAT  l/12X,»CHI  SQUARE  VALUE= • , F9.4 ) 

READ  (1,2)  NN,M 

N1=NN+1 

READ  (1,2)  1X1,1X2 

DO  33  1=1, M 

IY1=IX1*65539 

IF  (IY1)  51,61,61 

51  IYI=IY1+2147483647+1 

61  YFL1=IY1 
YFLl=YFLl*.4656613E-9 
IX1=IY1 
IY2=IX2*65539 

IF  (IY2)  52,62,62 

52  IY2=IY2+2147483647+1 

62  YFL2=IY2 
YFL2=YFL2*. 465661 3 E-9 
IX2=IY2 

Ri=YFLl 
R2=YFL2 

Z(1,I)=( (-2.*CL0G(R1) )**.5)*0C0S(R2*6.2832) 
33  Z(2,I)=( (-2.*CLCG(R1) )**. 5 ) *CS IN ( R2*6. 28  32 ) 
CM=M 
BN=NN 
B(l)=-3.5 
DO  10  1=2, Nl 

10  B(I  )  =  B(I-1)+0.14 
DO  9  K=l,2 

DO  9  1=1, NN 

9  C(K,I )=0.0 
DO  11  K=l,2 
DO  11  J=1,M 
H=-3.36 

DO  15  1=1, NN 
. IF  (Z(K,J).LE.M)  GO  TO  11 
15  H=H*0.14 

11  C(K,I )=C(K,I )41.0 
DO  14  1=1, Nl 

14  D(  I  )  =  B(T)/1.41A 
DO  18  K=l,2 
DD(K)=C.O 
AM(K)=C.O 
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SM.(K)-C.O 
'  DO    18    1=1, NN 

G(I  )=C.5*(DERF(D(  1  +  1)  J-DERFIDII)  D/CH 

FR(K,  l)=C(K,  !)/(CM*.14> 

DD(K)=CD<  KM  200*1  (FRIK,  I)-Gll))**2)/Gm 

AM(K)=C(K,l)*(B(I M0.071/CM+AMIK) 
18  SM(K)=SN(K)+(C(K,I )*( (B(I ) +0. 07 ) **2  )  /CM ) 

DO  30  K=l,2 

IF  (AN(K).EQ.O.O)  CO  TO  31 

DM(K)=SK(K)-(AM(K)**2) 

GO  TO  30 
31  DM(K)=SN(K) 
30  CONTINUE 

WRITE  (3,6) 

WRITE  (3,8)  (B(I),G(I),FR(1,I),FR(2,1),I=1,NN) 

DO  60  1=1,2 
WRITE  (3,42)  I 
WRITE  (3,40)  AtUDfCMCI] 
60  WRITE  (3,78)  CD ( I ) 
STOP 
END 
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C.3  PROGRAM  FOR  GENERATING  NEK  RANDOM  NUMBERS  (FUNCTION  OF  NORMAL 
RANDOf  PAIRS) 

IMPLICIT  INTEG£R**(I-N'J  t  REALMS  {  A-H,  C-Z  ) 

DIMENSION  C( 101 ) t B( ICl ItFRt 101 )tZl( 100OO) f Z2(IC0C0J 

2  F0RMA1  (2110) 

3  FORMAT  (2F7.2,3F5.2,F7.4) 

5  FORMAT  {/5(2X,«DX«,5X, 'FR',AX) J 

6  FORMAT  (1H1.1CX, • Vl= • , F7.2 , IX, « V2= • ,F7.2 , IX, «RHO= • ,F5. 2, IX, • Wl  = « , 
IF5.2,IX,*W2=« ,F5.2,1X,'ZTA=',F7.4) 

7  FORMAT  (/5( 1X,F6. 3, IX, F7. A, IX ) ) 

8  FORMAT  (10F0.5) 

40  F0RMAT(/12X,«MEAN=',F9.6t3X,«-  V  A  R  I  AN  CE  -  >',F9.6) 
•  REAO  (1,2)  NN,M 
READ  (1,2)  1X1,1X2 
DO  33  1=1, M 
IYI=IX1*65539 
IF(TYl)  51,61,61 

51  IY1=IYI+2147483647+1 

61  YFL1=IYI 
YFLl=YFLl*.4656613E-9 
1X1  =  1  Y  1 
IY2=IX2*65539 

IF  (IY2)  52,62,62 

52  IY2=IY2+2147483647+1 

62  YFL2=IY2 
YFL2=YFL2*. 46566 13E-9 
IX2=IY2 

Ri=YFLl 
R2=YFL2 

Zl(l)=( (-2.*DL0G(R1) )**.5)*DC0S(R2*6.28  32) 
33    Z2(I )=(l-2.*DL0G(Rl) ) **. 5 ) *DS INI R2*6.2832 ) 
CM  =  M 
BN=NN 
B(l)=-1.0 
DO    9    1=2, NN 

9  B(I  )  =  B(  I-D  +  0.02 
LS=l 

74    READ    (1,3)    VI, V2, RHC, Wl, W2, ZTA 

CF=DCCS(ZTA) 

SF=DSIN(ZTA) 

C(l)=C.O 

DO    10    1=2, NN 
10    C(I)=C.O 

SM=O.C 

AM=O.C 

DO    11    J=1,M 

Xl=Vl*Zl( JUKI 

X2=V2MRH0*Z1(J)+Z2(J)*DSCRT(1.-RH0*RH0) )+W2 

TT=DSCRT( l.+Xl*Xl+X2*X2) 

TR=(CF+X2*SF)/TT 
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H=-1.C 

DO  15  1=1, NN 

IF  (TR.LE.H)  GO  TO  il 

H=H+0.C2 
15  CONTINUE 
11  C(I  )=C(  D  +  1.0 

DO  30  1=1 |NN 

FR( I )=C(I »/(CP*.02J 

AN=(C(I)*B(I)/CM)+AN 

30  SH=(C(I)*(B(  I)**2  )/CM)+SM 
IF  (AN.EQ.G. J  GO  TO  3* 
0K=SP-(AM**2) 

GO  TO  32 

31  DK=SK 

32  WRITE  (3,6)  Vl,V2,RH0,Wi,W2,ZTA 
WRITE  (3,5) 

WRITE  (3,7)  (B( J) ,FR( J),J=i,NN) 

WRITE  (3,40)  AM, DM 

WRITE  (2,8)  (FR(J),J=l,NN) 

LS=LS+1 

IF  (LS.LE.6)  GO  TO  74 

STOP 

END 


Note:   1.  This  program  used  to  find  the  random  numbers  cos  61, 

2.  When  it  is  used  to  generate  random  numbers  of 
V 1  +  X2  +  Y2,  then  TR  =  TT,  B(l)  =  1.0,  B(l)  = 
B(I-1)  +  0.014.,  H  =  1.0,  H  =  H  +  0.0/j.,  FR(l)  = 
C(I)/(CM  «  .Ok). 

3.  When  it  is  used  to  generate  random  numbers  of 
l/Yl  +  X2  +  Y2,  then  TR  =  l/TT,  B(l)  =  0.0, 
B(I)  =  B(I-1)  +  0.01,  H  =  0.0,  H  =  H  +  0.01, 
FR(I)  =  C(I)/(CM  -::-  .01). 
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C.4  PROGRAM  FOR  THE  FIRST  KIND  OF  CURVE  FITTING 

IMPLICIT  1NTEGER*4(I-N)  ,REAL*8  «  A-H,  C-Z  ) 

DIMENSION  Y(51)tX(51),T(51,5)»PE(5i),Al6,6JtG(51),TPI5)fC(6), 
lB(5,5),TC(6),CRV(6),SV(5),D(2,6),CD(6),Fl5i),E(3),H{3,4),CA(3) 
DIMENSION  ED(IO) 
6  FORMAT  (4F8.4) 

8  FORMAT  (10X,2F9.5) 

9  FORMAT  (10X,  •  P= • , F7.3,6X,  •  E= • , F10.4 ,6X, ' KU=' ,  ID 

12  FORMAT  (10X, 'COEFFICIENTS') 

13  FORMAT  (IOX.8F12.5J 

14  FORMAT  (10X, «RTA=',F12.5) 
KSU=1 

15  READ,(Y( I) , 1=1,34) 
READ  (1,6)  XO,XL,P,DP 
NU=1 

MU=l 

P1=0.7C7 

N=4 

DX=(XL-X0)/33. 

X(l)=XC+DX/2. 

DO  1C  1=2, 34 

10  XII  )  =  XU-1)  +  DX 
DO  11  1=1,34 
T(I,l)=l.O 

DO  11  J=2,N 

11  T(I,J)  =  TU,J-1)*X(I  ) 
KT=1 

KU=1 
100  DO  30  1=1,34 
30  PE(I)sCEXP(-P*(X(I)-Pl)**2) 

00  40  1=1, N 

DO  40  J=1,N 

A(I,J)=0.0 

DO  40  K=l,34 
40  A(I,J)=A(I,J)4T(K,I )*T(K,J)*PE(K)*PE(K) 

MN=N+1 

DO  41  J=1,N 

DO  45  1=1,34 

45  C(I)  =  PE(I)*T(I,  J) 
DO  42  1=1,34 
TP(J)=C.O 

42  TP(J)  =  TP(J)+GU  ) 

CALL  INTG(G,34,DX,APA) 
41. AIMN, J)=ARA 

SY=O.C 

DO  46  1=1,34 

46  SY=SY+Y(I) 
DO  47  1=1, N 

47  A(IffMN)*SY*A(MNf I) 
A{MN,MN)=0.0 
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50 


55 


53 
54 


60 


65 


66 


68 


67 

85 

80 


70 


88 


89 


91 


300 


81 


DO  5C  1=1, N 

C(I)=C.O 

DO  50  J=l,34 

C(I  )=CU)+Y(  J)*T(  J,  [)*PE(J) 

C(MN)=1.0 

DO  55  1  =  1, N 

DO  55  J=l,N 

B(I,J}=A(I,J) 

IF  (N.EO.l )  GC  TO  53 

CALL  IiNVERS  (B,N) 

GO  TC  54 

BU,J)  =  1./B(1,1) 

DO  60  1  =  1, N 

TC( I)=C.O 

DO  60  J=1,N 

rcti  >  =  rcu  )+B(i,j)*cm 

CRF=0.0 

DO  65  1=1, N 

CRF=CRF+TC(l)*TP( I) 

DO  66  1=1, N 

CRV(  I)=CRF*A(MN,  l)*0.3 

CRV(Mf\)=0.0 

DO  68  1  =  1, N 

SV( I)=A(MN,I) 

DO  67  1=1, MN 

TC(i)=cm 

C(l )=C( IHCRVd  ) 

CALL  INVERS  (A,MN) 

KK=1 

DO  70  1*1, MN 

D1KK, I)=0.0 

DO  70  J=1,KN 

0(KK,I)=D(KK,I)+AU,J)*C(J) 

CALL  CCNST  ( MN, N, SV ,TP, 0,KK, TC,C ) 

KK=KK+1 

IFIKK.LT.3)  GC  TO  8C 

DO  8  8  1=1, MM 

DDU  )  =  CABS( (0(KK-1, I)/D(KK-2,I)  l-l.O) 

CRT=CC(1) 

DO  89  1=2, MN 

IF  ICRT.GE.ODU  )  )  GC  TO  89 

CRT  =CC( I) 

CONTIME 

IF  (CRT.LT.0.C001 )  GC  TO  300 

DO  91  1=1, MN 

D(1,I)=U(2, I) 

KK  =  2 

GO  TO  80 

DO  81  1=1,34 

FlI)=C.O 

DO  81  J=l,N 

F(I  )  =  F( I )+T(I,J)*0(2,J)*PE(I ) 

E(KT)=C.O 
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97 


899 

105 
400 
650 

600 
750 


620 


111 


DO  97 

E(KT>= 

WRITE 

WRITE 

WRITE 

WRITE 

WRITE 

ED(MU) 

HU=MU+ 

IF  (fl 

IF  (NO 

CALL  L 

MU=1 

KT  =  1 

NU=NU* 

GO  TO 

IF  (KL 

KT=KT+ 

IF  (KT 

P=P+DP 

GO  TO 

CCRT=( 

IF  (CC 

E(1)=E 

E(2)=E 

KT  =  3 

GO  TO 

T2=E(1 

IFIT2. 

KU=KU+ 

KT=3 

DO  605 

H(Itl) 

AI=l 

H(I,2) 

H(I,3) 

CALL  I 

DO  62C 

CA(I)= 

DO  62C 

CA(I)= 

P=-CA( 

GO  TO 

KSU=K$ 

IF  (KS 

STOP 

END 


1=1,34 

E(KT)+DABS(Y( I)-F(I) I 

(3,9)  P,E(KT),KU 
(3»12) 

(0(2,1), 1=1, N) 
0(2, KN) 
(Y( I ),F(I ), 1=1, 34) 


(3,13) 
(3,14) 
(3,8) 
=E(KT) 

1 

.LT.  6)  GO  TO 

•EQ.  3)  GO  TO 
SQP  (ED,P,DP) 


899 
111 


1 

100 

.EQ.2)  GO  TO  111 

1 

.GE.4)  GO  TO  400 

100 

E(KT-2)-E(KT-3) )* ( E ( KT-1 )-E ( KT-2 ) ) 

RT.LT.O.O)  GO  TO  600 

(2) 

(3) 

105 
)-2.*E(2)+E(3) 

LT.O.O)  GO  TO  650 
1 

1  =  1,3 

=  1.0 

=P-(3.0-Al )*0P 
=  H(I,2)*HU,2) 

NVERS  (H,3) 

1  =  1,3 
CO 

J=l,3 
H( I,J)*E( J)+CA( I ) 
2)/(2.*CA(3)  ) 
100 
0+1 
U.EQ.l)  CO  TO  15 
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SUBR01TINE  INTG(G,N,H,ARA) 
IMPLICIT  INTEGER*4U-N),REAL*8(A-H,C-Z) 
DIMENSION  G(34) 
ARA=G(1»-G(34) 
Nl=N-l 

DO  1C  1=2, Nl, 2 
10  ARA=ARA  +  't.G*GU  )+2.C*G(I+l) 
ARA=H*ARA/3.0 
!  RETURN 
END 

SUBROUTINE  CONST ( MN, N, SV,TP, 0, KK,TC,C ) 
IMPLICIT  INTEGER*4(I-N),REAL*8(A-HtC-Z) 
DIMENSION  SV(5),TPl5),D<2,6),TC(6)tCRV(6)tC<6) 
CRF=0.C 
DO  1C  1=1, N 
10  CRF=CRF+D(KK, I)*TP(I ) 
CRF=C(KK,MN)*CRF 
DO  20  1=1, N 
20  CRVd  )  =  CRF*SV(I  ) 
•   CRV(MN)=0.0 

DO  30  1=1, MN 
30  CII)=TC(I )+CRV( I ) 
RETURN 
END 

SUBROUTINE  LSCP  (ED.P.CP) 
IMPLICIT  INTECER*4(I-N)fREAL*8<A-H,C-Z} 
DIMENSION  ED(10),PT{10),C(5),B{3,3),A{3),0(3) 
PT(5)=P 
DO  1C  1=1,4 
AIM 
10  PT(I)=P-(5.-AI)*0P 
DO  20  1=1,5 
J  =  I-i 

cm=c.o 

DO  20  K  =  l,5 
20  C(I )=C(I)+PT(K)**J 

DO  30  1=1,3 

J=I-1 

D(I )=C.O 

DO  30  K=l,5 
30  D(I )=C(I )+ED(K)*(PT{K}**J) 

DO  40  1=1,3 

DO  4C  J=l,3 

K=J+I-1 
40  B(ItJ)=C(K) 

CALL  INVERS  (8,3) 

DO  50  1=1,3 

A(I)=C.O 

DO  5C  J=l,3 
50  All )=A( I )+B( I,J)*D( J) 

P=-A(2)/(2.*A(3) ) 

RETURN 

END 
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SUBROUTINE    [MVERS(A,NJ 

IMPLICIT    INTEGER*MI-N> , REAL* 8 ( A-H,C-Z ) 

DIMENSION  A(9,9),G(9),F(9) 

NN=N-1 

A(l,l)=1.0/A(l,l) 

00  21C  K=i,NN 

K  =  M+1 
250  DO  26C  1  =  1, M 

G(I)=C.O 

DO  26C  J  =  1,M 
260  G(I  )=G(I  )+A(I,J)*AUtK) 

Q=0.G 

DO  27C  1=1, M 
270  Q=Q+A(K,I)*Gll) 

Q=-Q+MK,K) 

A(K,K)-~1.0/Q 

DO  ?8C  1=1, M 
280  AU,K)=-G(I)*A(K,K) 

DO  29C  J=1,M 

F(JJ=C.O 

DO  29C  1=1, M 
290  F(J)  =  F(J)+AIK,I  J*A( I, J) 

DO  200  J=i,M 
200  A(K,J)=-F(J)*A(K',K) 

DO  21C  1  =  1, M 

DO  21C  J=1,M 
210  All, J)=A(I,J)-G(1)*A(K,J) 
•  RETURN 

END 
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C.5  PROGRAM  FOR -THE  SECOND  KINO  OF  CURVE  FITTINC 

DIMENSION  A<55,55),X(55)tYt2t55J,U(55),Vl55)»W(55)tXl(55,55)tR(55J 
1,T(55),DT(55),DM(55),DS(5),TY(55),C(55),8(55,55),DATA(2,98) 
DIMENSION  F(98) ,DF(98) 
860  FORMAT  ( lH0f5X, •COEFFICIENT* ■ /) 

865  F0RMAT(1H0,5X,' APPROXIMATE  VALUES^/) 

866  FORMAT  ( 1H0. 5X, • U' , / ) 

867  FORMAT  ( 1H0, 5X, • V • , / ) 
668  FORMAT  ( lHOf 5X, • N= ' , 14 ) 

869  F0RMAT11H  ,5X,«R'  ) 

870  FORMAT  (5E18.8) 

871  FORMAT  (IH1) 
879  FORMAT  (4E18.8) 

5  FORMAT  (10F8.5) 
7  FORMAT  (3F8.5) 
KP=1 
6C0  READ, NX, MK 

READ  (1,7)  (F(I), 1=1, 3) 
READ  (1,5)  (F(I ), 1=4, MX) 

00  9  1=1, MX 

IF  (F(I).LE.  0.0)  CO  TO  11 

GO  TO  9 
11  F(I)=C.0000i 
9  CONTINUE 

DO  10  1=1, MX 

DF( I )=0.4*SQRT(F( I ) ) 

DATA(1,I )=F( I)+DF(I ) 
10  DATA12, 1  )  =  F(I)-DFU  ) 

DO  666  KKA=1,MK 

READ,N 

CALL  SET(M,B) 

CALL  VECTOR  (DA TA,MX, M, Y ) 

Ml=M+l 

AM=M1 

S  =  0.0 

1  =  1 

15  U( I)=(Y(l,I )+Y(2,I) )/2.0 
V(I)  =  (Y(1,  I  )-Y(2, I ) )/2.0 
S=S*V(I)**2 
1=1  +  1 

IF(I.LE.M)  GO  TO  15 
S=SQRT(S) 

CALL  ARAMGE(B,U,M,A) 
KRITE  (3,871) 
N=M 

ND=N/2 
KL=-i 
NT=ND+2 
100  N=N*KL*ND 

IF(N.EC.M)GO  TO  8C0 
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25 


33 

30 


300 
38 


36 
37 


200 
40 

41 


45 


80 
70 


WRITE  (3,868)  N 

1  =  1 

TY(I)=C.O 

1=1  +  1 

IF(I.LE.M)  GO  TO  25 

1  =  1 

J=l 

Xlt I,J)  =  M I, J) 

J=J  +  1 

IF(J.LE.N)  GO  TO  30 

1  =  1  +  1 

I'FC  I.LE.M)  GO  TO  33 

KT=1 

KK=1 

KS=1 

1  =  1 

W(t)=l(I)+TY(I) 

1  =  1  +  1 

IF(I.LE.K)  GO  TO  33 

1  =  1 

C(I  )  =  C.O 

J=l 

C(I  )=C(  IJ+XKJ,  I  }*V)  i  J) 

J=J+1 

IF(J.LE.M)  GO  TO  37 

C(I)=C(I)*2.0/AM 

1=1  +  1 

IFU.LE.N)  GO  TO  36 

1  =  1 

R(I)=-U(I) 

J=l 

R(I)=R(I)+X1(I,J)*C(J) 

J=J  +  1 

IFU.LE.N)  GO  TO  41 

1  =  1  +  1 

IF(I.LE.M)  GO  TO  40 

WRITE  (3,869) 

SR=0.C 

1  =  1 

SR=SR+R( I)**2 

1  =  1  +  1 

IFU.LE.M)  GO  TO  45 

SR=SCRT(SR) 

IF(SR.LE.S)  GC  TO  70 

NT  =  N 

KL=1 

ND=NC/2 

IF(ND.GE.l)  GC  TO  ICO 

IF(N.LT.M)  GO  TO  80 

GO  TC  500 

ND=1 

GO  TC  100 

1  =  1 


at* 


72 


71 


233 

85 


95 


84 
86 


105 


115 


700 


710 


IF( 

GO 

1  =  1 

IF( 

KL  = 

ND= 

IF( 

IFl 

ND  = 

IF( 

GO 

1  =  1 

T(I 

IF( 

DT( 

GO 

DT( 

TY( 

GO 

TY( 

1  =  1 

IF( 

DS( 

1  =  1 

DS( 

1=1 

IFl 

KR1 

IF( 

KT= 

IFl 

GO 

DM( 

DMl 

KS= 

KT  = 

IFl 

GO 

R2  = 

Rl= 

R3= 

IFl 

DMl 

DM  I 

KS= 

IFl 

KK= 

GO 

KL= 

NF  = 

ND= 

IFl 

IFl 


ABSIRI  I  I  J.LE.VU  )  )  GO  TO  71 

TO  233 

+  1 

I.LE.M)  GO  TO  72 

-1 

ND/2 

ND.GE.l)  GC  TC  ICO 

N.GT.NT)  GO  TO  5C0 

1 

N.GE.2)  GO  TO  100 

TC  500 

)=RU)/VU) 

ABSITII ) J.LE.1.0)  GO  TO  95 

I  )=IABS(TU)  )-1.0)*RU) 

TO  84 

I)=C.O 

I  )  =  RU) 

TO  86 

l)  =  RU)-DTU  )*1.10 

+  1 

I.LE.M)  GO  TO  85 

KT)=0.0 


KT) 

+  1 

I.L 

TE 

OS  I 

KT  + 

KT. 

TO 

KS) 

KS) 

KS  + 

1 

KS. 

TO 

DMl 

DMl 

ABS 

R3. 

1)  = 

2)  = 
3 

KK. 
KK  + 
TO 

I 

N 

NC/ 

MD. 

N.t 


=  DS(KT)+ABSIDTU)  } 

E.M)  GO  TO  105 

13,870)  DS(KT) 

KTJ.LE.O. 1E-20)  GO  TO  555 

1 

GE.5)  GO  TO  115 

300 

=  DS(  1)+CS12)+CS13)*DS14) 

=DM{KS)/4.0 

1 

CE.4)  GO  TO  7C0 

3C0 

3) /DM12) 

2)/DM(l) 

IR1/R2-1.0) 

IT. 0.01)  GO  TC  710 

CM12) 

CM13) 

GT.M/A)  GO  TO  710 
1 

3C0 


GE.l)  GO  TO  ICO 
T.M)  GO  TO  72C 
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GO  TO  500 
720  ND=l 

GO  TO  100 
BOO  1=1 
BIO  C(I)=C.O 

J  =  l 
320  C(I)=C(I)+A( I,JJ*U(J| 

J=J  +  1 

IF(J.LE.M)  GO  TO  820 

Cm«C(I)*2.0/AM 

1  =  1  +  1 

IF(  I.LE.M)  GO  TO  810 

1  =  1 
830  TY( 1 )=C.0 

J=l 

840  TYII J=IY(I)+A(I, J)*C(J) 
J  =  J  +  l 

IF(J.LE.M)  GO  TO  84  0 
1=1  +  1 

IF(I.LE.M)  GO  TO  830 
GO  TO  765 
500  WRITE(3,868J  N 
1  =  1 

760  TY(I)=C.O 
J  =  l 

761  TY(I)=TYU)+X1(I,J)*C(J) 
J=J  +  1 

IF(J.LE.N)  GO  TO  761 

1  =  1  +  1 

IF(I.LE.M)  GO  TO  760 
765  WRITE(3,860) 

WRITE(3,870)  (C(I),I=l,N) 

WRITE(3,865) 
•  WRITE(3,870)  (TY ( I ) , 1=1 , M ) 

WRITE  (3,866) 

WRITE(3,870)  (U(I),I=1,M) 

WRITE  (3,867) 

WRITE  (3,870)  I Vf I  1 , I  =  1VM! 
666  CONTIME 

KP=KP4l 

IF  (KF.LE.2)  GO  TO  6C0 

STOP 

END 
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SUBROUTINE  VECTOR  t  H  ,  f-'l  , M,  Y  ) 

D I  MENS  I  ON  W I  2, 98 1 , Y (2 , 98 ) f  X ( 98 )  ,  Z ( 98  J 

IF(M.EC.Ml)  GC  TC  6C 

Al=Kl+l 

A2=M+1 

DO  1C  1=1, Mi 

AI  =  I 
10  X(I)  =  M/Al 

DO  15  1=1, M 

AI  =  I 
15  Z(l)  =  M/A2 

JM=1 

DO  30  1=1, M 
DO  35  J=JM,MM 
IF(Z(I).LT.X(J) )  GO  TO  35 
IF(Z(I).EQ.X(J) )  GC  TO  45 

CR=(Z(i)-x(j)  )*(zm-xu+i>) 

IF(CR.GT.O.O)  GO  TC  35 

Y(1,I)=W(1,J) 

Y(2,I)=W(2,J) 

GO  TO  65 
45  Y(l, I)=W(1, J) 

Y(2,I)=W(2,J) 

GO  TO  75 
35  CONTINUE 
65  JM=J 

GO  TO  30 
75  JM=J+l 
30  CONTINUE 

GO  TC  70 
60  DO  85  1=1, M 

DO  85  J=l,2 
85  Y(J, I)=W(J,I) 
70  RETURN 
•  END 
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SUBRGiriNC  ARANGE(A,U,K,C) 

DIMENSION  A (66, 66) ,U(66),B(66) ,KC(66) ,TP(66),C(66,66) 
200  FORMAT  1 1!U,5X,«  INKER  PRODUCTS'/) 
210  F0RMATU0E12.5) 
220  FORMAT  ( IHO, 5X, • CRCER  NO.'/) 

230  FORMAT (201 5) 

231  FORMAT  (1414! 
1  =  1 

10  B(I )  =  C.O 

J=l 
16  B(I  )=B(I)+A(J,I )*U(J) 

J=J  +  1 

IFU.LE.M)  GO  TO  16 

B(I)=ABS(8( I  )  ) 

1  =  1  +  1 

IF(  I.LE.M)  GO  TO  10 

G=-1.C 

WRITE  (3,200) 

L=l 

1  =  1 
20  J  =  l 

30  IF(BU).LT.G)  GO  TG  31 
IFIB(J).EC.G)  GO  TO  33 
G=B(J) 

L  =  J 

GO  TO  31 
33  IF(L.EC.l)  GO  TO  31 

IF(L.NE.KGU-l)  )  GO  TO  31 
L=J 

31  J=J  +  1 

IFU.LE.M)  GO  TO  30 

G=-l. 

KO(I)=L 

WRITE  (3,210)  B(L) 

B(L)=-2.0 

1=1  +  1 

IFU.LE.M)  GO  TO  20 

WRITE(3,??0) 

WRITE  (2,231)  (KC(I),  1  =  1, M) 

WRITE(3,230)  (KO (  I  ) , 1  =  1 ,M ) 

1  =  1 

35  K=KO(I) 
J=l 

36  C(J,  I)^A(J,K) 
J  =  J  +  1 

IFU.LE.M)  GO  TO  36 

1  =  1  +  1 

IFU.LE.M)  GO  TO  35 

RETURN 

END 
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SUOUCCTINE  SETCM,C1 

IMPLICIT  P.E  AL*8  ( A , 0 )  ,  REAL*4  1  B, Ct  E-H,CJ-Z )  ,  INT  EGER  I  I-  M  ) 

DIMENSION  C( ICO, ICC) 

Ml=M+l 

AM  =  M1 

AP=  3.  14109265358  979^+ (.32  38460-14) /AM 

Al=DSIN(AP) 

A2=CCCS(AP) 

ADl=Al 

AD2=A2 

CCltl)=SNGLlAl) 

1=2 
1C  AT=A2*AD1+A1*A02 

A2=-A1*A01+A2*A02 

A1  =  AT 

C(I»l)=SNCL(Al) 

AD3=AT 

A04=A2 

J=2 
20  AT=A2*AD3+A1*AD4 

A2=-Al*AD3+A2*AD4 

A1=AT 

C(I|J)=SNGL(A1) 

J=J  +  1 

IFtJ.LE.I )  GC  TO  20 

A1=AD3 

A2=AD'i 

1  =  1  +  1 

IF(I.LE.M)  GC  TO  10 

1  =  1 

40  11=1+1 
J=TI 

41  X(I  t J)-C(JV II 
J=J+l 

1FU.LE.M)  GO  TO  41 
1  =  1  +  1 

IF(I.LT.M)  GO  TO  40 
REIURN 
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A  study  of  probability  density  functions  of  the  angle,  6-^, 
between  the  surface  normal  and  the  direction  to  the  receiver  of 
a  rough  surface  in  a  radar  backscattering  problem  is  discussed. 
Both  the  analytical  and  the  Monte  Carlo  methods  for  finding  a 
probability  density  of  a  function  of  random  variables  are  intro- 
duced.  In  order  to  get  the  algebraic  expressions  of  the  density 
from  the  data  by  the  Monte  Carlo  method,  two  different  curve 
fitting  techniques  are  developed. 

Because  cos  6-,  involves  the  function  f  1  +  X2  +  Y2,  where 
X  and  Y  are  two  correlated  normal  random  variables,  the  problem 
of  finding  the  density  function  of  cos  Q-^   starts  from  finding 
the  densities  of  A/ \   +  X2  +  Y2  and  l/Vl  +  X2  +  Y2  with  a 
desire  to  find  the  former  density  by  the  information  method  con- 
tained in  the  latter  two.   The  Monte  Carlo  method  is  used  to 
find  the  densities  of  cos  9-,  with  each  different  parameter  CT , 
f,    and  6  due  to  the  failure  of  finding  the  density  function  of 
cos  6-.  analytically.   The  results  show  that  the  densities  of 
cos  6-i  are  approximately  Rice  distributions  centered  at  the 
value  of  cos  9^   when  CT   is  less  than  unity.   The  densities  of 
cos  0]_  are  also  discussed  by  changing  <J~ ,  f* ,    and  0. 


